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Abstract 

Effective  potential  energy  surfaces  (PESs)  are  calculated  for  the  nonadiabatic 
collision  5 (2P;-a)  +  H2(1'£g  ,v,j)  5(2Py^)  +  H2(1Zg,  v',/).  This  calculation 

employed  1  A',2  A  \  and  1  ~A "  adiabatic  PESs  numerically  calculated  at  the  state- 
averaged  multicon figurational  self-consistent  field  (SA-MCSCF)/configuration 
interaction  (Cl)  level  for  several  values  of  the  FE  bond  length,  FE  orientation  angle,  and 
boron  distance.  The  associated  nonadiabatic  coupling  terms  (NACTs)  were  calculated 
from  the  SA-MCSCF/CI  wave  functions  using  analytic  gradient  techniques.  A  line 
integral  through  the  NACTs  was  then  used  to  determine  the  adiabatic-to-diabatic  mixing 
angle  required  to  transform  from  the  1  A '  and  2  A '  adiabatic  basis  to  a  corresponding 
diabatic  basis.  When  all  nonadiabatic  coupling  terms  between  all  electronic  states  are 
considered,  the  line  integral  is  path  independent.  However,  only  NACTs  between  the 
1  A '  and  2  A '  states  were  considered  in  these  calculations,  and  the  line  integral  was 
therefore  path  dependent.  The  path  dependence  of  the  line  integral  was  used  to 
characterize  the  error  introduced  by  employing  a  truncated  set  of  adiabatic  states.  A 
method  for  reducing  the  effect  of  this  error  through  the  use  of  symmetry  derived 
boundary  conditions  was  developed.  The  resulting  diabatic  PESs  were  combined  with 
the  total  B  +  H2  rotational  kinetic  energy  and  boron  spin-orbit  coupling  to  yield  diabatic 
effective  PESs.  The  diabatic  effective  PESs  were  diagonalized  to  yield  adiabatic 
effective  PESs.  Diabatic  effective  PESs  data  was  extracted  for  the  equilibrium  Ho  bond 
length  and  used  to  calculate  inelastic  scattering  matrix  elements  using  the  time  dependent 


IV 


channel  packet  method.  These  matrix  elements  were  compared  previous  results  [D.  E. 
Weeks  et  al,  J.  Chem.  Phys.  125,  164301  (2006)]  to  observe  the  sensitivity  of  this 
calculation  to  the  input  electronic  structure  data. 
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THE  EFFECTIVE  POTENTIAL  ENERGY  SURFACES 


OF  THE  NONADIABATIC  COLLISION 


B(2PJa)  +  H2 (%+, v,;)  ~  S(2P;J  + 


I.  Introduction 


Motivation 

In  1926  Erwin  Schrodinger  published  his  landmark  paper  introducing  the  wave 
equation  which  would  later  bear  his  name.1  Schrodinger’s  equation  elegantly  reproduced 
the  results  achieved  by  Neils  Bohr  and  others  for  the  hydrogen-like  atom  and  was  quickly 
applied  to  other  atomic/molecular  systems.  However,  progress  on  solving  Schrodinger’s 
equation  for  systems  more  complex  than  the  hydrogen-like  atom  stalled.  Like  the  many- 
body  problem  in  classical  mechanics  the  resulting  equations  for  an  atomic/molecular 
system  could  not  be  solved  analytically. 

One  year  after  the  Schrodinger  equation  was  introduced,  Born  and  Oppenheimer 
proposed  that  the  dynamics  of  electrons  could  be  treated  separately  from  the  dynamics  of 
the  nuclei  in  atoms  and  molecules.  ’  This  approximation  was  motivated  by  the  fact  that 
electrons  move  much  faster  than  the  more  massive  nuclei.  This  approach  to  solving  the 
Schrodinger  equation  is  known  as  the  Bom-Oppenheimer  approximation  (BOA). 

Despite  the  utility  of  the  BOA,  in  1929,  Dirac  summarized  the  difficulties 
encountered  early  on  when  solving  the  Schrodinger  equation  for  atomic  and  molecular 


1 


systems  stating  the  following:  “The  underlying  physical  laws  necessary  for  the 
mathematical  theory  of  a  large  part  of physics  and  the  whole  of  chemistry  are  thus 
completely  known,  and  the  difficulty  lies  only  in  the  fact  that  the  exact  application  of 
these  laws  leads  to  equations  much  too  complicated  to  be  soluble.'’’’  However,  since  the 
advent  of  computers,  the  BOA  has  become  a  standard  tool  for  numerically  solving  the 
Schrodinger  equation  for  atomic  and  molecular  systems. 

In  the  BOA  the  positions  of  the  nuclei  are  fixed  by  setting  their  kinetic  energy 
operators  to  zero  in  the  Hamiltonian.  The  electron  kinetic  energy  operators  are  retained. 
The  resulting  electronic  Schrodinger  equation  is  solved  numerically  yielding  a  complete 
Hilbert  space  of  electronic  eigenfunctions,  and  corresponding  electronic  energy 
eigenvalues.  The  positions  of  the  nuclei  are  then  changed  and  a  new  electronic  problem 
is  solved.  This  process  is  repeated  for  a  variety  of  nuclear  configurations  using  nuclear 
coordinates  as  parameters  in  the  electronic  Hamiltonian.  The  resulting  collection  of 
ground  state  electronic  energy  eigenvalues  as  a  function  of  nuclear  coordinates  is  called 
the  ground  state  potential  energy  surface  (PES).  The  numerical  approaches  to  solving  for 
these  energies  and  wave  functions  include  Hartree-Fock  self-consistent  field  (HFSCF), 
configuration  interaction  (Cl),  and  multiconfigurational  self-consistent  field  methods.  In 
general  all  of  these  techniques  are  referred  to  as  ab  initio  methods. 

After  the  electronic  PESs  and  wave  functions  are  calculated,  the  nuclear  kinetic 
energy  operators  are  re-introduced  to  the  Hamiltonian.  The  resulting  nuclear/electronic 
Schrodinger  equation  is  solved  by  first  expanding  the  complete  nuclear/electronic  wave 
functions  in  tenns  of  the  electronic  wave  functions.  The  coefficients  of  this  expansion 
are  solved  by  representing  the  nuclear/electronic  Schrodinger  equation  in  the  electronic 
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basis.  This  is  called  the  adiabatic  representation.  In  this  representation  the  Schrodinger 
equation  has  terms  which  couple  electronic  eigenstates.  The  coupling  terms  are  referred 
to  as  nonadiabatic  coupling  terms  (NACTs).  The  BOA  ignores  the  NACTs  creating  a  set 
of  uncoupled  equations  equal  to  the  number  of  electronic  wave  functions  used  to  expand 
the  nuclear/electronic  wave  functions.4  Under  these  conditions  the  nuclear  dynamics  for 
a  given  electronic  eigenstate  are  influenced  by  a  single  PES  and  the  system  is  said  to 
behave  adiabatically. 

The  BOA  ignores  the  NACTs  whether  they  are  significant  or  not.4  The  NACTs 
can  be  large  when  the  electronic  eigenstates  have  a  strong  dependence  on  nuclear 
coordinates.  When  these  tenns  are  significant  the  adiabatic  representation  produces  a  set 
of  coupled  equations.  Under  these  conditions  the  nuclear  dynamics  for  a  given  electronic 
eigenstate  depend  on  multiple  PESs  and  the  system  is  said  to  behave  nonadiabatically. 

Nonadiabatic  systems  are  common,  especially  for  polyatomic  systems.5  Over  the 
past  several  decades,  a  great  deal  of  research  has  been  done  to  understand  nonadiabatic 
systems.6  Nonadiabatic  effects  play  an  important  role  in  the  radiationless  relaxation  of 
excited  electronic  states  and  the  photodissociation  of  molecules.7  These  effects  are 
recognized  to  play  a  critical  role  photobiochemical  processed  such  as  in  photosynthesis  in 
plants,  vision,  and  the  photochemistry  of  DNA.4  The  ability  to  accurately  model 
nonadiabatic  molecular  dynamics  is  applicable  to  many  areas  of  active  Air  Force 
research.  These  areas  include  the  modeling  of  collisional  transitions  between  atoms  and 
molecules  in  diode-pumped  alkali  lasers  (DPALs)  '  ,  and  the  modeling  high  energy 
density  materials  (HEDM)  for  explosives  and  fuel.11 
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Statement  of  Objectives 


Progress  in  modeling  nonadiabatic  systems  has  been  slow  due  to  the 
computational  requirements  they  impose.  Using  quantum  mechanical  formulations  '  . 
current  state-of-the-art  calculations  can  only  treat  the  nonadiabatic  dynamics  of  systems 
with  three  to  four  atoms.  Beyond  this,  classical  and  semi-classical  approximations19'25 
are  employed  to  achieve  results. 

This  research  uses  a  quantum  mechanical  treatment  to  explore  the  nonadiabatic 
effects  encountered  during  the  inelastic  collision  between  a  boron  atom  and  a  hydrogen 
molecule.  The  equation  for  this  scattering  event  is  given  by 

®(2p;J  +  (%W,/) «  b(2pj;)  +  h2(‘e+,v',/)  ( 1 ) 

In  this  interaction  the  boron  starts  in  a  “ P  electronic  state  and  is  free  to  transition  to  a 
different  P  state.  The  hydrogen  molecule  is  restricted  to  its  ground  electronic  state; 
however,  it  is  free  to  make  rotational  and  vibrational  transitions. 

The  interaction  of  atomic  boron  and  molecular  hydrogen  has  been  the  subject  of 
many  studies.  ’  As  stated  in  Weeks  (2006)  ,  general  interest  in  the  interaction  of 
atomic  boron  with  molecular  hydrogen  stems  in  part  from  the  high  degree  of  collisional 
quenching  that  is  observed  when  excited  metal  atoms,  including  boron,  react  with 
molecular  hydrogen.  Collisional  studies  of  excited  B(4 p  P )  with  H2  and  D2  have  also 

identified  a  significant  isotope  effect  on  reaction  rates  and  indicate  that  nonadiabatic 
mechanisms  may  influence  the  experimentally  observed  distribution  of  molecular 
products.  ’  Additional  interest  in  the  B+IB  system  is  derived  from  the  B+tB  van  der 

28  34  35 

Waals  complex  observed  in  molecular  beam  experiments'  ’  ’  and  in  boron  doped 

1131 

matrices  of  frozen  fB.  ’  Recent  interest  in  the  B  +  fR  system  stems  from  efforts  to 
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model  the  nonadiabatic  molecular  dynamics  of  the  inelastic  collision  B  (.2pi/2)  + 
H2(1T.g,j )  <-»  B(2P3/2 )  +  H2(1T.g,j' )  which  this  work  most  directly  extends.36'37  The 
“P  manifold  of  the  B  +  Tb  system  also  serves  as  a  proxy  for  the  excited  states  of  DPALs. 

The  research  presented  in  this  dissertation  has  five  main  aims.  First,  all  of  the 
theory  required  to  treat  the  interaction  given  by  Eq.  (  1  )  is  presented.  In  this  section  the 
theory  behind  the  BOA  is  presented  and  applied  to  the  B  +  H2  system.  The  method  for 
solving  the  Schrodinger  equation  when  the  BOA  breaks  down  is  presented.  The  specific 
case  of  an  atom  with  a  single  valence  electron  in  a  p  orbital  interacting  with  a 
homo  nuclear  diatomic  molecule  is  treated.  This  leads  to  the  method  for  calculating 
effective  potential  energy  surfaces  for  the  B  +  H2  system.  Finally  the  scattering  theory 
required  for  the  calculation  of  scattering  matrix  elements  using  the  Channel  Packet 
Method  (CPM)40'44  is  presented. 

Second,  ab  initio  calculations  of  the  three  adiabatic  PESs  corresponding  to  the 
boron  P  electronic  states  and  associated  NACTs  are  presented.  The  adiabatic  PESs  for 
this  work  were  calculated  at  the  state-averaged  multiconfigurational  self-consistent  field 
(SA-MCSCF)/configuration  interaction  (Cl)  level  by  Dr.  David  Yarkony  at  Johns 
Hopkins  Ehhversity.29, 45, 46  The  NACTs  were  calculated  from  the  SA-MCSCF/CI  wave 
functions  using  analytic  gradient  techniques.29, 47  Adiabatic  PESs  and  NACTs  were 
calculated  on  a  three  dimensional  grid  which  included  several  values  of  the  molecular 
hydrogen  bond  length.  Previous  ab  initio  calculations  of  the  adiabatic  PESs  for  the  B  + 
H2  system  fixed  the  molecular  hydrogen  bond  length  at  the  equilibrium  value  of  1 .402 
a.u.  ’  ’  ’  This  restriction  forces  the  hydrogen  molecule  to  remain  in  its  ground 
vibrational  state.  The  data  presented  in  this  dissertation  are  a  prerequisite  for  future  work 
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in  which  vibrational  transitions  of  the  hydrogen  molecule  are  included  in  the  calculation. 
Furthennore,  previous  work  did  not  directly  calculate  the  NACTs.'  ’  The  nonadiabatic 
effects  were  modeled  using  other  methods.  In  the  past  the  direct  calculation  of  NACTs 
was  avoided  primarily  because  NACTs  can  become  singular  for  some  systems.  These 

40 

singularities  make  it  very  difficult  to  solve  the  Schrodinger  equation  numerically. 
However,  methods  which  use  the  NACTs  are  able  to  make  an  estimate  of  overall 

48 

accuracy  of  the  results  that  other  methods  cannot  make. 

Third,  the  adiabatic  representation  is  transfonned  via  a  unitary  transfonnation  to 
form  the  diabatic  representation.  When  this  transformation  is  chosen  properly  the 
NACTs  are  negligible  in  the  diabatic  representation,  eliminating  any  singularities  in  the 
NACTs.  The  form  of  the  Schrodinger  equation  is  also  greatly  simplified  in  the  diabatic 
representation.  This  work  employs  a  line  integral  through  the  vector  field  defined  by  the 
NACTs  to  arrive  at  the  adiabatic-to-diabatic  mixing  angle  used  in  the  unitary 
transformation.49’ 50  The  symmetry  of  the  B  +  H2  system  allows  the  value  of  this  mixing 
angle  to  be  predicted  for  certain  configurations.  ’  This  allows  the  accuracy  of  the 
diabatic  PESs  to  be  estimated.  This  work  also  presents  a  new  method  for  reducing  the 
source  of  error  introduced  by  restricting  the  electronic  basis  to  include  only  the  boron  P 
electronic  states. 

Fourth,  the  diabatic  PESs  are  cast  into  a  form  best  suited  for  wave  packet 
propagation.  This  is  done  by  fitting  the  diabatic  PESs  to  reduced  Wigner  rotational 
matrix  elements  and  employing  an  interpolation  scheme  to  yield  fit  coefficients  at  any 
desired  nuclear  coordinate.'  ’  The  fit  coefficients  are  used  to  construct  effective  PESs 
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which  include  the  effects  of  spin-orbit  coupling  and  rotational  dynamics.  This  work 
presents  both  diabatic  and  adiabatic  forms  of  these  effective  PESs. 

Finally,  the  data  corresponding  to  the  equilibrium  bond  length  of  the  hydrogen 
molecule  is  extracted  from  these  effective  PESs  and  the  CPM  is  used  to  compute 
scattering  matrix  elements.  These  scattering  matrix  elements  are  compared  with  previous 
one  dimensional  scattering  results.  ’  ’ 
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II.  Theory 


The  Hamiltonian 

In  the  framework  of  quantum  mechanics,  the  Hamiltonian  operator  is  used  to 
detennine  the  time-evolution  of  the  state  of  a  system  as  part  of  the  time-dependent 
Schrodinger  equation  (TDSE).  The  general  form  of  the  TDSE  is  given  by 

ik^  =  H\'V)  (2) 

at 

where  IT1)  is  a  wave  function.  While  the  wave  function  IT1)  can  be  chosen  arbitrarily  it  is 
often  chosen  as  a  superposition  of  eigenstates  of  the  Hamiltonian  determined  by  the  time- 
independent  Schrodinger  equation  (TISE).  The  general  form  of  the  TISE  is  given  by 

H\<P)  =  E\<P)  (3) 

where  |<b)  is  an  eigenstate  of  the  Hamiltonian  operator  H.  E  is  the  energy  eigenvalue 

associated  with  the  eigenstate  |<t>).  The  set  of  all  eigenstates  of  the  Hamiltonian  fonn  a 
complete  set  of  orthogonal  states  which  can  be  used  to  express  other  states. 

The  Hamiltonian  of  a  system  of  interacting  atomic  nuclei  and  electrons  is  given 


by 


H  -  fN  +  fe  +  VNN  +  VeN  +  Vee  (  4  ) 

The  operator  fN  is  the  sum  of  all  the  kinetic  energy  operators  for  the  nuclei. 


•  -ft 

x 

a= 1 


Pa 

2mn 


(5) 


In  Eq.  (  5  )  the  sum  is  carried  out  over  the  index  a  up  to  nn,  the  total  number  of  nuclei  in 
the  system.  The  mass  of  each  nuclei  is  given  by  ma.  The  momentum  operator  for  each 
nuclei  is  given  by  pa. 
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The  second  operator  in  Eq.  (  4  ),  fe,  is  the  sum  of  all  the  kinetic  energy  operators 


for  the  electrons. 


'-I!; 

i=  l 

In  Eq.  (  6  )  the  sum  is  carried  out  over  the  index  i  up  to  the  number  of  electrons  in  the 
system  given  by  ne .  The  mass  of  each  electron  is  the  same  and  is  given  by  me .  Like  the 
nuclei,  each  electron  has  its  own  momentum  operator  which  is  given  by  pi- 

The  third  operator  in  Eq.  (  4  ),  VNN,  is  the  sum  of  the  Coulomb  potential 
interaction  operators  between  nuclei. 


ZaZpe2 


\  1  \  1  fi 


a=l  (3>a  1  “ 

In  Eq.  (  7  )  the  sum  is  over  all  unique  parings  of  nuclei  in  the  system.  The  charge  of  an 
electron  is  e  and  the  charge  numbers  of  each  nuclei  are  given  by  Za  and  Zo.  The  position 
operators  of  each  nuclei  are  given  by  qa  and  qp.  The  operator  q  represents  the  set  of 
position  operators  {qx,  q2, ... },  one  for  each  coordinate  used. 

The  fourth  operator  in  Eq.  (  4  ),  VeN,  is  the  sum  of  the  Coulomb  potential 
interaction  operators  between  nuclei  and  electrons. 


=  y  y  zae2 

r—i  Z—t  I  qi  —  qc 


In  Eq.  (  8  )  the  sum  is  carried  out  over  each  nuclei  and  electron  pair,  where  a  is  used  to 
index  nuclei  and  i  is  used  to  index  electrons. 

The  last  operator  in  Eq.  (  4  ),  Vee,  is  the  sum  of  the  Coulomb  potential  interaction 
operators  between  electrons. 
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ne  n . 


(9) 


V„ 


q  i-qj 


In  Eq.  (  9  )  the  sum  is  carried  out  over  all  unique  pairings  of  electrons  in  the  system.  The 
operators  in  Eqs.  (  5  )  through  (  9  )  are  given  in  abstract  form.  The  specific  form  of  the 
operators  will  depend  on  the  representation  chosen.  When  these  equations  are  combined, 
the  Hamiltonian  H  for  a  system  of  nn  nuclei  and  ne  electrons  can  be  written  (in  atomic 
units)  as51 


Tip  Tin 


a= 1  1=1  a=lp>a'^a  ^P' 

ne  ne 


£  =  1  a= 1 


Qi-  Qi 
1=1  1>1  n 1 


(10) 


In  atomic  units  the  mass  and  charge  of  an  electron  are  e  =  1  and  me  =  1.  The  Plank 
constant  h  is  also  unity  in  atomic  units.  The  number  of  terms  in  Eq.  (  10  )  increases 
rapidly  as  electrons  and  nuclei  are  added  to  the  system.  The  B  +  H2  system  has  three 
nuclei  and  seven  electrons.  The  Hamiltonian  describing  this  system  will  have  55  tenns. 
The  TISE  for  the  B  +  Hz  system  cannot  be  solved  analytically.  Numerical  methods  for 
solving  the  TISE  in  this  form  for  systems  more  complex  than  a  single  hydrogen  atom  are 
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still  being  developed. 


The  Bom-Oppenheimer  Approximation 

In  1927,  Bom  and  Oppenheimer  published  their  approach  for  solving  the  TISE  for 
a  molecule.  Their  approach  was  motivated  by  comparing  the  relative  masses  of 
electrons  and  protons.  Electrons,  being  roughly  1800  times  less  massive  than  a  proton, 
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typically  move  much  faster  than  the  atomic  nuclei.  Thus,  the  electrons  can  be  thought  of 
as  moving  through  the  electric  field  produced  by  stationary  nuclei.  The  nuclei,  in  turn, 
can  then  be  thought  of  as  moving  through  the  average  electric  field  produced  by  the 
quickly  moving  electrons.  This  observation  motivated  Bom  and  Oppenheimer  to  treat 
the  dynamics  of  the  electrons  separately  from  the  nuclear  dynamics. 

Under  the  assumption  that  the  nuclei  are  stationary,  the  nuclear  kinetic  energy 
operators  are  eliminated  from  Eq.  (  10  )  defining  a  new  Hamiltonian: 


tip  n„  n- 


H, 


elec 


=y^+yy  ZaZ*  +yy  Ja 

ii2  tk.p?a\qa~qp\  9al 


(ii) 


Eq.  (  1 1  )  is  the  electronic  Hamiltonian.  As  in  Eq.  (  10  )  a  and  /?  index  nuclei  while  i  and 
j  index  electrons.  The  nuclear  positions  are  now  parameters  in  Helec  rather  than 
operators.  The  symbol  q  represents  the  set  of  coordinates  {qlt  q2, ... }.  The  set  of  nuclear 
and  electronic  coordinates  are  represented  by  qn  and  qe  respectively.  The  solutions  to 
the  electronic  TISE  have  a  parametric  dependence  on  nuclear  positions.  The  electronic 
TISE  is  given  by 


Helec\j(<ln))  =  KleM n)  I K<ln))  (  12  ) 

where  \j(qn))  refers  to  the  jth  eigenstate  of  the  electronic  Hamiltonian  with  the  nuclear 

positions  for  a  given  configurations  of  nuclei  given  by  the  set  of  nuclear  coordinates  q. 

The  eigenstate  \jXqn))  is  expressed  in  the  electronic  coordinate  representation  as 

(qe\j((ln))  =  (13) 

In  Eq.  (  13  )  the  dependence  of  <t>y  (qn,  qe )  on  nuclear  and  electronic  coordinates  is 

shown  explicity. 
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In  Eq.  (  12  )  £'g;ec(qn)  refers  to  the  /h  energy  eigenvalue  of  the  corresponding 
eigenstate.  In  the  electronic  Hamiltonian  the  nuclear  repulsion  tenn  VNN  becomes  a 
constant  for  a  given  nuclear  configuration  and  is  added  to  the  energy  eigenvalue 
Fgjec(q,n)  to  yield  the /h  PES.  The  PES  for  a  given  eigenstate  is  a  function  of  nuclear 
coordinates.  The  dimensionality  of  the  PES  increases  as  the  number  nuclear  degrees-of- 
freedom  (DOF)  increases. 

For  each  nuclear  configuration  there  are  an  infinite  number  of  electronic 
eigenstates.  These  states  form  a  set  of  orthonormal  functions  which  is  complete,  as 
shown  in  Eq.  (  14  ). 


<*(<!n)l/0?n)> 


J  dqe(i(qn)\qe)(qe\j(qn)) 
I  dqe^t(qn,qe)<t>j(qn,qe) 


Sa 


(14) 


The  eigenstates  for  different  nuclear  configurations  are  not  orthogonal  to  one  another  in 
general. 

Even  after  formulating  the  electronic  Hamiltonian  Heiec,  the  corresponding 
electronic  TISE  must  be  calculated  numerically  using  methods  such  as  the  Hartree-Fock 
(H-F)  approximation,  configuration  interaction  (Cl),  and  multi-configuration  self- 
consistent  field  (MCSCF)  methods.  These  methods  are  referred  to  as  ab  initio  techniques 
(Latin:  from  the  beginning  ).  There  are  many  software  packages  available  that  perform 
these  calculations.  Commonly  used  programs  include  GUASSIAN,  GAMESS, 

MOLPRO,  and  COLUMBUS.  These  programs  numerically  solve  the  electronic  TISE  by 
creating  wave  functions  for  each  electron  (spatial  orbitals)  and  combining  them  in  linear 
combinations  to  fonn  molecular  orbitals.  The  linear  combination  of  spatial  orbitals  is 
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optimized  using  variational  methods.51  The  resulting  wave  functions  are  the  eigenstates 
of  Helec- 

Using  Eq.  (  1 1  )  the  full  nuclear  and  electronic  USE,  Eq.  (  3  ),  can  be  rewritten  as 

£m  =  (tN  +  ttetejm  =  E\V)  (15) 

The  set  of  all  electronic  eigenstates  form  a  complete  set  and  can  be  used  as  a  basis  to 

represent  other  functions.  The  solutions  to  the  full  TISE  can  be  represented  as 

W<?„)>  =  2  Fj(qa)  ly(fln)>  (  16  ) 

j 

In  Eq.  (  16  )  the  full  nuclear/electronic  wave  function  is  expanded  using  electronic 
eigenstates  where  Fj(jqn )  is  the  expansion  coefficient  and  is  a  function  of  nuclear 
coordinates.  The  coefficients  Fj(qn)  can  be  solved  by  substituting  Eq.  (  16  )  into 
Eq.  (  15  ),  multiplying  from  the  left  by  (i(qn)  |,  and  integrating  over  electronic 
coordinates. 


(i(Qn)\{TN  +  Helec  -E)  I^F;(qn)  |y(<fn)>] 


(17) 


f  i  ) 

This  representation  is  called  the  adiabatic  representation.  When  electronic  eigenstates  are 
used  to  represent  the  solutions  of  the  full  TISE  they  are  referred  to  as  adiabatic  states  and 
their  corresponding  PESs  are  called  adiabatic  PESs. 

Since  the  total  energy  eigenvalue  if  is  a  constant,  Eq.  (  17  )  can  be  further 
simplified  by  using  the  orthonormal  property  of  the  electronic  eigenstates  (  Eq.  (  14  )  ). 


{i(qn)\(-E)\^Fj(qn)  \j(qn))  J 


-E^FjiqJSy 


-EFiCqJ 


(18) 
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The  electronic  Hamiltonian  Helec  does  not  contain  derivatives  with  respect  to  nuclear 
coordinates,  thus  it  has  no  effect  on  Fj(qn).  This  results  in  the  term  derived  in  Eq.  (  19  ). 


elec )  X  Fj(q n)  ly(<7n)>  |  =  F}elec(qn)  ^  T;(qn)5j 


(19) 


=  £,iiec(9n)Ti(qn) 


The  operator  T N  involves  derivatives  with  respect  to  nuclear  coordinates,  and  since  both 
the  coefficient  F)(qn)  and  the  eigenstate  17 (<7n))  depend  on  nuclear  coordinates,  the 
product  rule  must  be  used  when  evaluating  this  term.  For  clarity,  the  coordinate 
representation  is  used  to  express  this  inner  product. 


(Kqn)\{TN  )  |^V;(qn)  \j(qn)) 

=  -U  -j-  [  dqe(pf*(qn,qe)V^n{Fj(qn)(pf(qn,qe )} 
z  /  i  Trla  j 
J,cc 

=  _2  AirT J  d‘3'echf*(9n<9e){fKqn)Vqn<h/(qn<9,e)  +  2  (v^Fy (<?n)) 


(20) 


],a 

•  mfle))  +  ^/(fliv  <7e)V^nFy  (<?„)} 

The  superscript  4  has  been  added  to  the  coordinate  representation  of  the  eigenstates 
e)  to  clearly  identify  them  as  adiabatic  states.  The  nuclear  coordinate 

dependence  of  the  del  and  Laplacian  operators  V£?n  and  are  indicated  by  the  subscript 
qn.  The  form  of  Eq.  (  20  )  can  be  simplified  by  defining  the  following  terms: 


»5  =  j 

(21) 

ka~- 
lJ  2. 

[  dqe<pf*(qn,qe)V2qn<pf(qn,qe) 

(22) 
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The  terms  f  fj  and  fcfj  couple  electronic  eigenstates  via  their  dependence  on  nuclear 


coordinates  and  are  called  derivative  and  kinetic  coupling  terms,  respectively.  They  are 
referred  to  as  nonadiabatic  coupling  terms  (NACTs).  The  derivative  coupling  terms 
(DCTs)  are  vector  quantities  while  the  kinetic  coupling  terms  are  scalar  values.  The 
DCTs  have  a  vector  component  for  each  nuclear  coordinate  used  to  describe  the  system. 
When  referring  to  a  specific  component  of  the  vector  field  defined  by  the  DCTs  the 
notation  is  used,  where  q  denotes  the  specific  nuclear  coordinate.  In  the  literature 


the  following  shorthand  notation  for  DCTs  is  frequently  encountered: 

=  {rfl&t) 

The  DCTs  are  antihennitian  as  expressed  in  Eq.  (  24  )  (see  Appendix  A). 

tA  —  —rA* 

Lij  Lji 


(23) 


(24) 


Combining  the  results  of  Eqs.  (  1 8  )  —  (  20  )  and  using  the  definitions  in  Eqs.  (  21  )  and 
(  22  ),  Eq.  (  17  )  can  be  expressed  as 


+  2 ■  V,„  +  2 KfiFjiqJ 


(25) 


+  {Elelec  ~  E)Fi(q n)  =  0 


Eq.  (  25  )  now  represents  the  TISE  in  the  adiabatic  representation.  Aside  from  treating 
the  electronic  problem  separately  from  the  complete  nuclear/electronic  problem,  no  other 
approximations  have  been  made.  However,  expanding  the  TISE  using  Eq.  (  16  )  has 
created  a  system  of  coupled  differential  equations.  When  the  NACTs  are  either  small 
enough  to  be  negligible  or  simply  neglected  the  TISE  becomes  uncoupled  and  can  be 
written  as 


(-^a^  +  EleieM)Fi(qn)  =EFi(qn ) 


(26) 
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This  is  referred  to  as  the  adiabatic  Born-Oppenheimer  approximation.6  In  Eq.  (  26  )  the 
expansion  coefficient  Fj(qfn)  is  interpreted  as  the  nuclear  wave  function  whose  dynamics 
are  governed  by  a  single  PES  FeZecfan)-  The  tenn  Eleiec(qn )  is  called  an  adiabatic  PES. 

In  later  calculations  VA  fan)  E elec  fan)  will  be  used  as  an  alternate  notation  for  the 

adiabatic  PES  associated  with  the  zth  electronic  eigenstate  |i(qrn)). 

If  a  system  is  accurately  modeled  by  Eq.  (  26  )  the  system  is  said  to  behave 
adiabatically.  However,  not  all  systems  can  be  modeled  accurately  by  the  adiabatic 
BOA.  When  the  NACTs  are  too  large  to  be  neglected  the  adiabatic  BOA  is  invalid.  The 
generalized  Hellmann-Feynman  theorem,  Eq.  (  27  )  ,  gives  the  clearest  indication  of 
when  a  system  will  behave  nonadiabatically. 

->A  _  {'■(cIn)\VqnHelec |)fan))  t  27  1 

lJ~  P/ fan) -P,"  fan) 

This  theorem  is  derived  in  Appendix  B.  When  ( i(qn )  |Vqn  Heiec\j(qn))  is  non-zero  and 
the  difference  between  VjA(qn)  and  PP*fan)  is  small  f fj  will  be  large.  In  the  case  where 
P/fan)  =  Vt\q  n),  the  derivative  coupling  f fj  will  be  singular.  Points  where  two  PESs 
intersect  and  rf  is  singular  are  called  conical  intersections.5  Systems  that  have  conical 
intersections,  like  the  B  +  H2  system,  cannot  be  modeled  using  Eq.  (  26  ). 

The  most  common  approach  to  solving  Eq.  (  25  )  with  non-negligible  coupling 
tenns  involves  transfonning  from  the  adiabatic  basis  to  a  new  basis,  the  diabatic  basis,  in 
which  the  coupling  terms  are  negligible.49’ 50  However,  practical  computational 
limitations  restrict  the  number  of  electronic  eigenstates  that  can  be  considered.  To  reduce 
the  computational  overhead,  ab  initio  data  is  calculated  only  for  a  restricted  subspace  of 
the  full  electronic  Hilbert  space  consisting  of  the  electronic  states  most  strongly  coupled 
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S'} 

to  one  another.  Restricting  the  subspace  neglects  any  weak  coupling  between  states 
within  the  subspace  and  external  states  and  can  limit  the  effectiveness  of  the  adiabatic-to- 
diabatic  transformation  (ADT). 54-56  These  considerations  can  be  best  illustrated  by 
examining  the  B  +  H2  system. 

The  B  +  H?  System 

2  2 

In  its  ground  state,  boron  has  a  ls~2s"2p  electronic  state.  The  unpaired  2 p 
electron  gives  rise  to  a  P  spectroscopic  term  consisting  of  three  degenerate  atomic 
orbitals  energetically  separated  from  other  boron  orbitals.  As  the  boron  atom  approaches 
the  hydrogen  molecule  the  degeneracy  is  lifted,  however,  they  remain  energetically  close 
to  one  another.  By  Eq.  (  27  )  it  is  assumed  that  these  orbitals  will  be  strongly  coupled  to 
one  another  while  only  weakly  coupled  to  states  outside  the  P  manifold.  Consequently, 
the  electronic  basis  is  truncated  to  include  only  the  three  states  within  the  ~P  manifold. 

In  general,  the  B  +  FE  system  has  Cs  =  {1,  <7}  symmetry  where  1  is  the  unit 
operator  and  0  corresponds  to  reflection  though  the  B,  H2  plane.  The  symmetry  of  the 
orbitals  is  defined  as  follows: 

A'  =>  cr|<T>)  =  l|<b> 

A"  =>  u|0)  =  — 1|<J>)  (28) 

According  to  Eq.  (  28  )  orbitals  that  are  unaffected  by  a  reflection  about  the  B,  FE 
plane  have  A '  Cs  symmetry.  Figure  1  shows  an  example  of  this  type  of  orbital. 
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Figure  1.  An  example  of  a  boron p  orbital  with  A'  Cs  symmetry 


In  Figure  1,  the  green  circles  represent  the  hydrogen  atoms,  the  blue  circle  represents  the 
boron  atom,  the  blue  lobe  represents  the  portion  of  the  p  orbital  with  positive  phase,  and 
the  red  lobe  represents  the  portion  of  the  p  orbital  with  negative  phase. 

If  a  reflection  about  the  B,  FF  plane  flips  the  orientation  of  the  orbital  then  the 
orbital  has  A  "  Cs  symmetry.  An  example  of  this  type  of  orbital  is  shown  in  Figure  2. 


Figure  2.  An  example  of  a  boron p  orbital  with  A  "  Cs  symmetry 


In  Figure  2  a  reflection  about  the  B,  H2  plane  would  cause  the  lobes  of  positive  and 
negative  phase  of  a  p  orbital  to  change  positions.  The  two  boron  2 p  orbitals  with  A '  Cs 
symmetry  are  labeled  as  1  A'  and  2  A '  with  1  ~A "  labeling  the  orbital  with  A "  Cs 
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symmetry.  These  orbitals  have  corresponding  adiabatic  PESs  V1 2 Ar,  V2  2A>,  and  VA«. 
Here  the  superscript  A  has  been  omitted  to  avoid  confusion  with  the  symmetry  label. 
Given  that  the  B  +  H2  involves  a  single  atom  interacting  with  a  homonuclear  diatomic 
molecule  the  adiabatic  PESs  satisfy  the  symmetry  relation  V(r,n  —  9,  R)  =  V (r,  9,  R ). 

The  TISE,  Eq.  (  25  ),  for  the  B  +  H2  system  with  a  truncated  electronic  basis  can 
be  written  as 


u-\y<2 

2  Z_ 1  m„ 


IV 


2  /—>  m„ 


-S-/ 

2  L—i  171  a  J 


YJ^2-yqn  +  <]  o\ 


J 


(29) 


(V**  0 

+  0  V22A, 

Vo  0 


(E  0  0\/fW\ 

0  E  0  1 1  f  2  24'  J 

Vo  0  eJ  V  fa»  J 


In  Eq.  (  29  )  the  NACTs  have  been  separated  from  the  kinetic  energy  operators  and  PESs. 
Orbitals  with  different  C5  symmetry  will  not  couple  giving  rise  to  the  off-diagonal  zeros 
in  the  NACT  matrix.  This  is  discussed  in  Appendix  C.  Although  the  eigenstate  with  A  " 
Cs  symmetry  does  not  couple  with  the  other  eigenstates  it  is  still  included  in  the 
calculation.  All  three  functions  span  the  boron  2 p  subspace  and  are  required  when 
transfonning  from  a  basis  labeled  by  the  angular  momentum  quantum  numbers  to  one 
labeled  by  Cartesian  coordinates.  As  a  result  of  Cs  symmetry,  the  matrix  containing  the 
coupling  tenns  is  block  diagonal  with  a  2-by-2  block  and  a  1  -by- 1  block. 

The  coupling  matrix  can  be  simplified  further  by  choosing  real  valued  eigenstates. 
Given  the  antihennitian  property  of  the  DCTs,  real  valued  functions  forces  the  diagonal 
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coupling  terms  to  be  zero.  This  in  turn  forces  the  diagonal  kinetic  coupling  terms  to  be 
zero  due  to  the  following  relationship  between  derivative  and  kinetic  coupling  tenns. 


y  '  fA  _  A  _  =>A 
VQn  LIJ  ~  l]  LIJ 


4j  =  V  •  +  ifj  ■  rfj 


(30) 


According  to  Eq.  (  30  ),  if  the  derivative  coupling  tenn  for  a  given  pair  of  eigenstates  is 
exactly  equal  to  zero  everywhere  the  corresponding  kinetic  coupling  term  will  also  be 
equal  to  zero. 


The  B  +  H?  Coordinate  System 

The  PESs  as  well  as  the  NACTs  depend  on  the  nuclear  coordinates  used  to 
describe  the  system.  The  internal  DOF  of  the  B  +  FE  system  are  described  by  Jacobi 
coordinates  R ,  the  distance  between  the  boron  atom  and  the  center-of-mass  (CM)  of  the 
hydrogen  molecule;  r,  the  bond  length  of  the  hydrogen  molecule;  and  6,  the  polar  angle 
fonned  by  the  axis  containing  the  boron  atom  and  the  axis  containing  the  hydrogen 
molecule  which  intersect  at  the  CM  of  the  hydrogen  molecule.  Figure  3  illustrates  the 

nn 

same  body  fixed  (BF)  and  space  fixed  (SF)  coordinate  system  used  by  Weeks  et  al.  A 
body  fixed  (BF)  coordinate  system  is  chosen  in  which  the  boron  atom  and  the  FE  CM  are 
fixed  on  the  BF  z  axis.  The  orientation  of  the  hydrogen  molecule  is  determined  by  the 
angle  6  and  the  projection  of  the  FE  molecular  axis  on  the  BF  xy  plane.  The  azimuthal 
angle  0  is  formed  between  the  BF  x  axis  and  this  projection.  This  angle  is  not  shown  in 
Figure  3. 


20 


z 


Figure  3.  Body  fixed  and  space  fixed  coordinate  system  used  to  describe  the  B  +  H2  system37 

The  BF  axis  is  oriented  in  the  SF  coordinate  axes  by  using  Euler  angles  a  and  ft  and 
constraining  the  BF  y  axis  to  lie  in  the  SF  XY  plane.  The  Euler  angle  a  is  formed  by  the 
projection  of  the  BF  z  axis  onto  the  SF  XY  plane  and  the  SF  X  axis.  The  Euler  angle  /?  is 
the  polar  angle  formed  between  the  BF  z  axis  and  the  SF  Z  axis.  The  origin  of  both  the 
BF  and  SF  axes  coincide  at  the  B,  H2  CM.  In  Figure  3  the  CM  is  displaced  away  from 
the  boron  atom  for  clarity.  In  this  coordinate  system  the  1  ~A '  and  2  ~A '  orbitals  are 
chosen  to  correspond  to  the  pz  and  px  orbital  of  the  boron  atom  respectively,  while  the  A  " 
orbital  corresponds  to  the  py.  ’ 

The  Adiabatic-to-Diabatic  Transfonnation 

The  TISE  for  the  B  +  FE  system  in  the  adiabatic  representation  is  given  by  Eq. 

(  29  ).  The  NACTs  make  it  difficult  to  solve  Eq.  (  29  )  directly.  As  mentioned  earlier,  it 
is  desirable  to  transform  from  the  adiabatic  representation  to  the  diabatic  representation  in 
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which  the  NACTs  are  negligible.  Many  methods  for  calculating  diabatic  states  have  been 
described  in  the  literature.57  These  methods  can  be  divided  into  two  groups:  those  that 
calculate  the  adiabatic-to-diabatic  transformation  (ADT)  using  the  NACTs,  and  those  that 

58 

detennine  the  diabatic  states  using  the  matrix  elements  of  various  Hennitian  operators. 
The  methods  in  the  latter  group  avoid  the  computational  overhead  of  calculating  the 
NACTs  over  a  grid  in  nuclear  configuration  space  dense  enough  to  capture  the  strong 
variation  NACTs  can  experience  as  a  function  of  nuclear  coordinates.59'63  This  variation 
is  especially  severe  near  conical  intersections,  where  the  DCTs  become  singular.5 
However,  it  is  generally  accepted  that  methods  using  NACTs  are  more  accurate  in 
calculating  diabatic  states.48' 54 

With  the  truncated  B  +  H2  basis  used  in  Eq.  (  29  )  the  ADT  is  defined  as 
cos  Y(qn)  -sin  y(qn)  0 

sin  y(qn)  cos  y(qn)  0  1  NVy  (31) 

0  0  V  \  <V  ) 

where  the  angle  y(qn)  is  the  ADT  mixing  angle.57  The  adiabatic  states  are  labeled 
according  to  their  symmetry.  The  diabatic  states  are  given  the  labels  xx,  yy,  and  zz.  In 
matrix  fonn  Eq.  (  31  )  is  given  by  <bD  =  (Lfid>A.  The  matrix  U  represents  the  diabatic-to- 
adiabatic  transformation.  The  ADT  transfonnation  is  unitary,  thus  ILfilU  =  1.  As 
discussed  in  Appendix  C,  only  states  with  A '  Cs  symmetry  will  mix;  hence,  the 
transformation  matrix  HJt  has  the  same  block  diagonal  form  as  the  coupling  matrix  in  Eq. 

(  29  ).  The  ADT  mixing  angle  y(qn)  is  calculated  for  each  nuclear  configuration. 

The  derivative  coupling  tenn  between  1  A '  and  2  A '  states  in  the  adiabatic 
representation  can  be  transfonned  into  the  diabatic  representation  by  using  Eq.  (  3 1  )  to 
make  the  following  substitutions  into  Eq.  (  21  ). 
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<*>l2A'((ln.<le)  =  c*,zz0?n.<le)cOSyG?n)  +  ®xx(Rn>  Re)  sin YCRti) 

(32) 

^22A'(Rn>Re)  =  <bzz(qfn,qe)siny(qfn)  +  ^xx(qn.Re)  cosy(qn) 

After  simplifying  the  resulting  expression  Eq.  (  33  )  is  obtained.  The  details  of  this 
derivation  are  given  in  Appendix  D. 

T12  =  Vqny(<7?l)  +  t£,  (  33  ) 

Under  the  assumption  that  the  DCTs  in  the  diabatic  representation  are  equal  to  zero  Eq. 

(  33  )  simplifies  to 

Vqny(qn)  =  (  34  ) 

For  systems  with  a  single  nuclear  coordinate  this  equation  can  be  integrated  directly  to 
yield  the  ADT  mixing  angle.  For  systems  with  more  than  one  nuclear  coordinate  the 
single  integral  becomes  a  path  integral  through  the  vector  field  defined  by  each 
component  of  the  DCTs.50’ 64-66  For  the  case  of  B  +  FE  the  integral  takes  the  following 
form: 

r  e  r 

Y(R,B,r)  =  j  T$il2\rogodR  +  j  Tg12\R^d6  +  j  T^12\Rgdr  +  y0  (35) 

Ro  ro 

Eq.  (  35  )  assumes  a  specific  order  of  integration  for  the  Jacobi  coordinates.  Eq.  (  35  ) 
gives  the  ADT  mixing  angle  y  to  an  overall  constant  yQ.  A  reference  point  (r0, 90l  R0)  in 
the  Jacobi  coordinate  space  is  chosen  to  start  the  line  integral  for  all  detenninations  of  the 
ADT  mixing  angle.  The  ADT  mixing  angle  for  this  reference  point  is  set  to  zero  or  some 
arbitrary  constant. 

Under  the  assumption  that  ff2  =  0  the  value  of  the  mixing  angle  should  be  path 
independent  as  indicated  by  Eq.  (  36  ). 

Vqn  X  *12  =  0  (  36  ) 
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Eq.  (  36  )  follows  directly  from  the  vector  calculus  theorem  that  states  that  the  curl  of  the 
divergence  of  a  vector  field  must  be  zero:  Vqn  X  (vqnY(qn)')  =  0.  When  this  condition 


is  met  the  ADT  completely  eliminates  the  NACTs  and  the  resulting  diabatic  surfaces  are 
called  strictly  diabatic  surfaces.48’ 64, 67 

Under  the  assumption  that  the  ADT  successfully  eliminates  the  NACTs,  the  TISE 
given  by  Eq.  (  29  ),  becomes 

y -\y-  °  o  \ 
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The  resulting  diabatic  PESs  are  given  by  the  following  equations. 

Vxx(r,  9,  R)  =  Vy2 Af  (r ,  6,  R)  cos2  y(r,  6,  R)  +  V22Ai(r,  6,  R)  sin2  y(r,  6,  R) 
Vzz(r,  9,  R )  =  Vl2 y(r,  B,  R )  sin2  y(r,  6,  R)  +  V22 A>(r,  9,  R)  cos2  y(r,  9,  R) 
Vxz(r,9,R )  =  sin  y(r,  9,  R)  cos  y(r,  9,  R)  (VyA'(r,9,R)  -  V22 Ar(r,9,R)) 
Vyy(r,9,R)  =  VA’i(r,  9,  R) 


(37) 


(38) 


However,  as  a  consequence  of  truncating  the  electronic  basis  set  and  keeping  only  those 
terms  in  the  2P  manifold,  DCTs  involving  states  outside  the  basis  have  been  neglected. 
This  introduces  error  into  the  ADT  mixing  angle.  This  is  best  understood  by 
decomposing  the  derivative  coupling  vector  field  into  a  transverse  and  longitudinal 

/TO 

component. 


yA  _  yAX  yA.l 
l12  ~  l12  '  l12 


(39) 


The  transverse  and  longitudinal  component  have  the  following  properties. 


V  '  =  0 


(40) 


24 


(41) 


With  these  definitions,  Eq.  (  33  )  becomes 

V<Zn  X  *12  =  Vqn  x  *12 

The  error  associated  with  the  neglected  coupling  tenns  gives  rise  to  a  transverse 
component  of  the  derivative  coupling  vector  field  which  introduces  path  dependence  in 
the  line  integral  used  to  calculate  the  ADT  mixing  angle/  '  This  error  cannot  be 
removed  and  is  often  referred  to  as  nonremovable  component  of  the  derivative  coupling 
field.57  While  Eq.  (  34  )  can  still  be  solved  for  the  ADT  mixing  angle,  the  error  will 
accumulate  through  each  path.  Diabatic  surfaces  affected  by  this  error  are  tenned  quasi- 
diabatic  surfaces.  The  accuracy  of  the  resulting  quasidiabatic  states  is  determined  the 
magnitude  of  the  residual  coupling.  In  the  cases  where  this  residual  coupling  is  small  the 
truncated  basis  is  expected  to  adequately  capture  the  nonadiabatic  effects  in  the 
dynamics.69 

Near  a  conical  intersection  the  nonremovable  coupling  is  small  relative  to  the 
removable  coupling,  making  a  smaller  overall  contribution  to  the  line  integral.70  As  the 
distance  from  the  conical  intersection  increases  the  removable  and  nonremovable 
couplings  approach  the  same  magnitude.  Consequently,  the  path  used  for  the  line 
integral  has  typically  been  restricted  to  regions  of  nuclear  configuration  space  near  a 
conical  intersection  to  minimize  the  relative  contribution  of  nonremovable  error  to  the 
path  integral.72'76  The  restricted  range  of  nuclear  coordinates  also  minimizes  the  errors 
introduced  by  finite  grid  spacing,  thus  allowing  the  location  of  the  conical  intersection  to 
be  determined  with  a  high  degree  of  accuracy.  Xu  and  coworkers  have  extended  the 

line  integral  technique  to  search  for  conical  intersections  throughout  all  nuclear 
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configuration  space  for  H3.  However,  their  method  for  calculating  the  NACTs  relies  on  a 
semiempirical  determination  of  the  ADT  matrix  in  conjunction  with  Hellmann-Feynman 
theorem  and  ensures  a  curl-free  vector  field.77  Thus  their  line  integral  is  constrained  to  be 
path  independent. 

Another  approach  for  solving  for  the  ADT  mixing  angle  involves  taking  the 
divergence  of  Eq.  (  33  )  to  obtain 

%n  •  =  V^n7(qn)  +  v,n  •  fj>2  (  42  ) 

Using  the  definitions  in  Eq.  (  40  )  and  assuming  that  the  resulting  diabatic  NACTs  are 

negligible  Eq.  (  42  )  simplifies  to 

=  -V  '  4'  (43) 

Eq.  (  43  )  is  a  Poisson’s  equation  which  allows  the  ADT  mixing  angle  to  be  solved 

without  perfonning  a  line  integral.  Work  has  been  done  to  solve  Eq.  (  43  );  however,  the 

boundary  conditions  are  unknown  and  must  be  approximated.55-  78  Furthermore  there  are 

numerical  restrictions  on  the  regions  where  Eq.  (  43  )  can  be  solved. 

The  Asymptotic  Basis 

The  calculation  of  scattering  matrix  elements  is  not  conveniently  done  by  solving 
the  TISE  in  the  diabatic  electronic  representation  (  Eq.  (  37  )  ).  The  time-dependent  CPM 
is  one  approach  for  calculating  these  scattering  matrix  elements.43’44  In  this  method  the 
Hamiltonian  for  the  B  +  H2  system  is  represented  using  states  defined  in  the  asymptotic 
limit  of  the  system. 

In  the  asymptotic  limit  is  reached  as  the  distance  R  between  B  and  Hi  becomes 
large.  In  the  asymptotic  limit  the  B  +  H2  system  consists  of  two  non-interacting  systems: 
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the  boron  atom  and  the  H2  molecule.  The  boron  atom  has  a  nucleus  and  five  electrons 
and  the  tb  molecule  has  two  nuclei  and  two  electrons.  The  terms  in  the  full  Hamiltonian 
can  be  classified  according  to  the  system  they  belong  to.  The  full  Hamiltonian  for  the  B  + 
H2  system,  excluding  spin-orbit  coupling,  can  be  written  as 


H  =  f £  +  f '"2  +  TeB  +  re[  +  ff2  +  Vef  +  Vel  +  Eoff  (  44  ) 

The  tenns  TB  and  feB  represent  the  kinetic  energy  operators  for  the  nucleus  and  electrons, 


^  H  w  h 

respectively,  of  the  boron  system.  The  terms  TN  2  and  Te  2  represent  the  same  for  the  H2 
system.  The  term  represents  the  Coulomb  interaction  potential  between  particles  in 
the  boron  system  only.  The  term  Vel  2  represents  the  same  for  the  H2  system.  The  tenn 
Vet  represents  the  Coulomb  interaction  potential  between  the  boron  system  and  the  H2 
system.  This  term  has  the  form  (in  atomic  units) 


Vt 


el 


-2s 

hj  1 


7  7 


tmj  I 


(45) 


In  Eq.  (  45  ),  the  subscript  a  refers  to  particles  associated  with  the  atom,  while  the 
subscript  m  refers  to  particles  associated  with  the  molecule.  The  atomic  number  of  a 
particle  is  given  by  Z.  The  set  of  position  operators  for  each  charged  particle  is  given  by 
q.  Finally,  the  last  term  E0f  j  represents  a  constant  energy  offset  which  has  no  effect  on 
the  eigenstates. 

In  the  asymptotic  limit  of  large  R  the  contribution  of  the  electrostatic  interaction 
potential  to  the  total  energy  goes  to  zero,  and  the  asymptotic  Hamiltonian  becomes 


Ho  =  Tn+  t"2  +  HBei  +  +  Eoff  (  46  ) 

The  eigenbasis  of  the  asymptotic  Hamiltonian  H0  is  a  direct  product  of  the  states 

associated  with  each  term  of  Eq.  (  46  ).  The  term  HBt  is  the  electronic  Hamiltonian  of 
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the  boron  system.  This  operator  is  the  sum  of  all  the  Coulomb  interaction  tenns  between 
charged  particles  in  the  boron  system  as  well  as  all  the  electron  kinetic  energy  operators 
for  the  boron  system.  It  can  be  expressed  as 


H?l  =  ?eB  +  Vet  =y^  +  I|/A 
l  2me  fjlRai-Ctaj 


(47) 


Hel2  represents  the  same  for  the  H2  system  and  is  expressed  as 

Pm,i  \ '  ^mi^vij 


fjH 2  "pn2  .  t >n2  \  '  V m,i  \ 1 

"«■  ~T‘  +v“  -Z^AZ 


e  |  Qmi  Qmj  | 


(48) 


As  discussed  earlier,  the  basis  of  electronic  eigenfunctions  was  truncated  to  include  only 
the  2 p  orbitals  associated  with  the  unpaired  boron  electron.  These  orbitals  are  also 
eigenfunctions  of  the  term  HBt.  The  electronic  eigenstates  of  the  boron  atom  are 
represented  by  |T,  A).  When  spin-orbit  coupling  is  ignored  the  states  \  £  =  1,  A  =  0,  +1) 

are  degenerate  with  an  energy  of  E B  corresponding  to  the  boron  fP  term.  As  the  other 
sources  of  angular  momentum  of  the  B  +  H2  system  are  coupled  the  electronic  eigenstates 
for  the  boron  atom  will  be  indexed  by  spin-orbit  coupling  labels. 

The  asymptotic  basis  of  the  electronic  hydrogen  molecule  Hamiltonian  includes 
only  the  ground  electronic  state  represented  by  |  ).  In  the  asymptotic  limit,  its  energy 

eigenvalue  is  given  by  the  vibrational  potential  energy  curve  VH2  (r)  which  depends  on 
the  hydrogen  bond  length  r.  These  eigenstates  fonn  the  asymptotic  electronic  basis 
which  is  represented  by  |X,  A)  =  1 1Sg )  (x)  |T  =  1,A  =  0,±1). 

Representing  Eqs.  (  47  )  and  (  48  )  in  this  basis  yield  the  asymptotic  energy  of  the 
system,  as  shown  by 
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(S', +  H"i2} M  =  EbpS^Sa>a  +  VH2(r)S^S 


a'a 


(49) 


The  operators  //®  and  H^2  are  diagonal  in  this  representation.  The  potential  VHz  (r) 
represents  the  vibrational  potential  energy  of  H?  in  its  ground  electronic  state  |  ). 

In  previous  work  ’  ,  the  hydrogen  bond  length  was  restricted  to  its  equilibrium 
value  r  =  req.  When  r  is  held  at  a  fixed  value  the  asymptotic  energy  for  the  hydrogen 
molecule  becomes  a  constant  and  can  be  eliminated  along  with  the  asymptotic  energy  of 
the  boron  atom  by  choosing  E0ff  appropriately.  In  this  work  the  requirement  that 
r  —  req  is  relaxed.  By  relaxing  the  restriction  on  r,  the  asymptotic  energy  of  ffi  is  no 
longer  a  constant  with  respect  to  r.  Consequently  the  r  dependence  of  VH2  (r)  cannot  be 

2  p  m  _  _  2  p 

eliminated  and  the  constant  energy  offset  E0^  =  —EB  is  used  to  eliminate  EB  only. 

In  anticipation  of  representing  the  full  Hamiltonian  in  the  asymptotic  electronic 
basis  1 2,  A),  Eq.  (  46  )  can  be  written  as 

H0  =  Tg  +  f"2  +  H"2  (  50  ) 

Using  the  BF  and  SF  coordinates  developed  earlier,  the  nuclear  kinetic  energy  operators 

can  be  expressed  in  a  CM  frame  as  79' 80 


fB  ,  f  tf2  _  PR  ,  Pr  J2  L2 

1 N  '  1 N  —  ~  r  “  >"  -  —  -r 


(51) 


2Hb,h2  2[1h2  2  2hb,h2P2 

In  Eq.  (  51  )  the  reduced  masses  are  given  by  E-BHl  =  2mgm/;2/(mB  +  2 mH2)  and 

/ iH2  =  mH/ 2.  The  terms  pR  and  pr  are  momentum  operators  conjugate  to  the  coordinates 
R  and  r.  The  angular  momentum  operator  ]  corresponds  to  the  rigid  rotor  comprised  of 
the  H2  molecule.  The  angular  momentum  operator  L  corresponds  to  the  rigid  rotor 
comprised  of  the  boron  atom  and  the  CM  of  the  H2  molecule. 
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The  asymptotic  Hamiltonian  H0  can  now  be  expressed  as 


Pr  !  Pr  |  J2  |  £2 

2  Pb,h2  2  /%2  2/%2r2  2  nBH2R2 


+  Hel 


(52) 


The  operators  in  Eq.  (  52  )  can  be  grouped  into  a  vibrational,  translational,  and  rotational 
piece. 


Pr  ,oh2 

^r+H" 


(53) 


Eq.  (  53  )  describes  the  vibration  of  H2  in  its  ground  electronic  state.  This  piece  will  be 
expressed  in  the  coordinate  representation  using  the  Jacobi  coordinate  r. 

Eq.  (  54  )  describes  the  translation  of  the  boron  atom  with  respect  to  the  CM  of 


H^trns  _ 
0  — 


2  Pb,h2 


(54) 


The  eigenfunctions  of  this  translation  operator  are  plane  waves.  Each  plane  wave  is 
labeled  by  its  corresponding  momentum  (in  atomic  units)  given  by  |k). 

Eq.  (  55  )  describes  the  rotation  of  the  H2  molecule  and  the  tumbling  motion  of 
the  boron  atom  around  the  Ho  molecule. 


2  PhJ2  2Pb,h2r2 


(55) 


There  are  four  contributions  to  the  total  angular  momentum  J  for  the  B  +  H2  system:  j,  L, 
f,  and  5.  These  angular  momentum  must  be  coupled  to  produce  a  set  of  quantum 
numbers  to  index  the  rotational  eigenfunctions.  Dubernet  and  Hutson  examine  five 
different  coupling  schemes.79  This  work  uses  the  coupling  scheme  labeled  “case  1  A”  in 
their  paper.  This  same  scheme  is  described  by  Jouvet  and  Beswick.81  This  coupling 


scheme  is  illustrated  in  Figure  4. 
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Figure  4.  Angular  momentum  coupling  scheme  "case  1A"  described  by  Dubernet  and  Hudson74 

In  this  scheme  the  angular  momenta  £  and  s  associated  with  the  orbital  angular  moment 
and  spin  of  the  unpaired  boron  electron  couple  to  give  ja,  which  has  a  projection  o>  on  the 
BF  z  axis.  Earlier  in  this  section  the  electronic  eigenstates  of  boron,  neglecting  spin,  were 
indexed  by  their  orbital  angular  momentum  A.  In  the  coupled  angular  momentum 
scheme  they  are  indexed  by  the  spin-orbit  indices  ja  and  co.  The  angular  moment  of  the 
H2  rotor  j  has  a  projection  k  onto  the  BF  z  axis,  and  is  not  coupled  with  ja.  The  total 
angular  moment  is  given  by 

]  =  L  +  j  +  )a  (  56  ) 

The  total  angular  momentum  is  chosen  to  have  a  projection  on  the  BF  z  axis  given  by 

P  =  k  +  co  and  a  projection  M  onto  the  SF  z  axis.  In  the  centrifugal  sudden  (CS) 
approximation  P  —  k  +  <u  is  held  constant.  The  angular  momentum  L  of  the  tumbling 
motion  of  the  boron  atom  with  the  CM  of  the  hydrogen  molecule  has  no  projection  along 
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the  BF  z  axis  since  the  boron  atom  and  the  CM  of  the  hydrogen  molecule  lie  along  the  BF 
z  axis. 


Using  these  labels,  the  basis  functions  of  the  spin-orbit  coupled  BF  basis  are 
labeled  as  follows: 

'Hmi-Ms)  <57> 

The  label  /  is  defined  to  correspond  to  properly  chosen  values  for  the  quantum  numbers 
in  Eq  (  57  ).  Spin-orbit  coupling  of  the  unpaired  boron  atom  leads  tol/2  and  3/2  as  the 
possible  values  of  ja,  allowing  co  to  have  values  ranging  from  —3/2  to  3/2  in  integer 
steps.  There  are  an  infinite  number  of  rotation  energy  states  for  the  H2  rotor,  allowing/ 
to  range  from  zero  on  up  in  integer  steps.  The  projection  k  ranges  from  —j  to  j  in  integer 
steps.  There  is  no  quantum  number  associated  with  the  tumbling  motion  of  the  boron 
atom  and  the  CM  of  the  hydrogen  molecule  in  this  coupling  scheme.  The  possible  values 
of  the  total  angular  momentum  /  range  from  1/2  on  up  in  integer  steps. 

The  complete  basis  of  the  asymptotic  Hamiltonian  is  given  by  the  direct  product 
of  the  eigenfunctions  for  each  of  the  tenns  in  Eq.  (  52  ). 


K.r.R)  =  10®  k>®  |K>®  l1^) 
=  MPkts»'T‘R’*) 


(58) 


The  eigenfunctions  form  an  orthononnal  set  as  shown  by 


(C,r',R'\Z,r,R)  =  8?j8(r' -r)8(R' -  R) 


(59) 


Since  the  hydrogen  molecule  is  restricted  to  the  ground  state  it  will  be  eliminated  from 


the  notation.  These  functions  are  now  used  to  represent  the  full  Hamiltonian. 
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The  Asymptotic  Representation 

The  full  Hamiltonian,  this  time  with  a  spin-orbit  coupling  tenn  included,  will  be 
represented  in  the  asymptotic  given  by  Eq.  (  58  ). 

F 


-  Pr  Pr  J2 
H  =  +  -P—  +  -1 — T  +  - 


.  .  ~  1  _■  +  Hi2  +  Vei  +  H  so 

2-Pb,h2  2  Hh2  2  Hh2t  2 Pb,h2 ^ 

In  Eq.  (  60  ),  Hso  is  the  spin-orbit  term  and  the  constant  energy  offset  E0^  has  been  used 
to  eliminate  the  contribution  of  the  electronic  boron  Hamiltonian  H Jj.  Since  the 
asymptotic  basis  functions  are  not  eigenfunctions  of  the  full  Hamiltonian,  the  fonn  of  the 
TISE  in  this  representation  will  not  be  diagonal. 

There  are  an  infinite  number  of  states  in  the  asymptotic  basis.  For  practical 
computational  reasons  this  basis  must  also  be  truncated.  For  this  work,  only  the  states 
where  J  =  1/2,  M=  +  1/2,  and  P  =  +  1/2  are  included.  By  truncating  the  electronic 
basis  functions  to  include  only  the  2P  states  of  the  boron  atom  £  —  1,  s  —  1/2,  and 
ja  =  1/2,  3/2.  The  hydrogen  molecule  is  restricted  to  the  ground  electronic  state.  As 
shown  in  Eq.  (  61  ),  the  basis  indices  (  Eq.  (  58  )  )  can  be  simplified  by  eliminating  the 
labels  that  remain  constant. 


(60) 


K.r.R) 


J  J  a 
k  a) 


,r,R 


(61) 


The  matrix  elements  of  the  first  tenn  of  the  Hamiltonian  (  Eq.(  60  )  ),  the  radial  kinetic 
energy  term,  are  given  by  (in  atomic  units) 

-1  d2 


j  Ja  u' 

,  ,  r 

\k  CO 


Pr 


2  Pb,c 


J  Ja 


:R8v 


kco’r'R  2  nB,HRdR2“"xx 


(62) 


where  8X'X  is  given  by 

s,’,  =  -  ow  -  r') 


(63) 
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The  matrix  elements  of  the  vibrational  kinetic  energy  term  are  given  by  Eq.  (  64  ). 


7  Ja  o' 

,,  ,<  r  >  K 

\k  Cl) 


A 2 

Vr 
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7  ja  r  n\  _  ^  ^  r  s; 

k(o'  ’  2 jiHrdr2  xx 


(64) 


The  matrix  elements  of  the  H2  rotational  energy  are  given  by  Eq.  (  65  ). 


7  Ja  n' 

, ,  r  < K 

k  (x) 
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J  Ja 

k  CO 


,r,R 


7(7  +  1) 
2^h  r 


2  x  x 


(65) 


Eqs.  (  62  ),  (  64  ),  and  (  65  )  represent  tenns  of  the  Hamiltonian  that  are  diagonal  in  this 
basis. 


The  matrix  elements  of  the  angular  momentum  associated  with  the  tumbling 
motion  of  the  boron  atom  and  the  hydrogen  molecule  CM  is  not  diagonal  in  this 
representation.  The  matrix  elements  are  obtained  by  using  Eq.  (  56  )  to  express  L2  as 


L2  =  j2  +  j2  +  j2  -  2/  •  j  -  2]  •  ja  +  2]  •  ja 

=  J2  +J2  +jl  ~  J+J _  -J-J+  ~  2 jzJz  -J+Ja-  —  J~ja+  ~  2 JJaz  +  J+Ja- 
+  J-Ja+  +  VzJciy 


(66) 


In  Eq.  (  66  )  j+  and  J  are  the  raising  and  lowering  operators  of  the  total  angular 
momentum  operator.  Likewise,  ]+,}-,  Ja+  and  ja-  are  the  raising  and  lowing  operators 
of  the  angular  momentum  operators  j  and  ja  respectively.  The  matrix  elements  can  then 
be  calculated  using  the  established  properties  of  raising  and  lowering  operators  with  one 
exception.  The  operation  of  j-  are  reversed  as  a  consequence  of  the  BF  angular 
momentum  calculus.  ’  The  tenns  in  Eq.  (  66  )  can  be  classified  as  either  diagonal  or 
off-diagonal  terms.  The  terms  diagonal  in  this  representation  are  those  that  do  not 
include  raising  and  lowering  operators.  These  terms  operate  on  the  basis  functions  in  the 
following  way 
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(67) 


[f  +  f  +  Ja  -  2/Jz  -  2 ]zJaz  +  2jzJaz}  |  V,  R 

=  {/(/  +  1)  +j(j  +  1)  +ja(Ja  +  1)  -  2(/c  +  ca)2  +  2/cca} 
The  off-diagonal  terms  are  given  by  Eqs.  (  68  ),  (  69  ),  and  (  70  ). 
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The  matrix  elements  of  the  tumbling  kinetic  energy  tenns  are  then  given  by 
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(70) 


+  {(/(/  +  1)  -  P(P  +  D)(70  +  1)  -  k(k  +  1))}1/2  x  S^S^S^Sir  -  r')8(R  -  R ') 

+  {(70  +  1)  -  P(P  +  D)(7a0a  +  1)  -  k(k  +  D)}1/2  x  SJ^MS^kSJ},^S(r  ~  r'WR  ~  R ') 
+  {(70  +  1)  “  ^(fe  +  l))(7a0a  +  1)  “  &>("  +  1))}  7 

X  -  r  w  -  *')] 


(71) 


Eq.  (  71  )  is  grouped  into  four  terms,  each  of  which  couples  states  for  which  Aj  —  0  and 


Aja  —  0.  The  second  and  third  terms  couple  states  for  which  AP  =£  0.  It  is  common  to 
neglect  these  terms  by  invoking  the  CS  approximation.  "  This  approximation  assumes 
that  the  interaction  takes  place  rapidly  enough  that  the  projection  P  of  the  total  angular 
momentum  is  conserved.  This  work  uses  the  CS  approximation. 
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The  spin-orbit  coupling  operator  is  given  by  Hso  —  <f?  •  s  where  ?  and  s  are  the 
orbital  and  spin  angular  momentum  operators  of  the  unpaired  2pl  boron  electron.  The 
matrix  elements  of  this  operator  in  the  asymptotic  basis  are  given  by  Eq.  (  72  ). 

{k'u"r'’R'\Rso\JkJy>R)  =  \{ja{ja  +  1}  -  ^  -  <s  +  !)}«*'*  (  72  ) 

In  this  expression,  the  atomic  boron  value  of  <f  is  4.876  x  10  5  atomic  units.  For  this 
calculation  the  spin-orbit  coefficient  <f  is  assumed  to  be  constant  and  is  set  to  the  atomic 
boron  value.  This  is  called  the  “pure  precession”  approximation.  In  general,  ^  changes 
as  a  function  of  the  nuclear  coordinates  due  to  the  perturbation  of  the  boron  atomic 
orbital. 

The  electrostatic  interaction  potential  VeL  has  a  straightforward  representation  in 
the  coupled  angular  momentum  asymptotic  basis  when  first  expanded  using  a  standard 
multipole  expansion.90  The  multipole  expansion  has  the  fonn  (in  atomic  units) 

.j  _  \  ’  ZaiZmj 
Bl  ^  |7ai  “  Rmj\ 

^  '  (73) 

In  Eq.  (  73  )  the  electrostatic  interaction  potential  is  expanded  in  terms  of  renormalized 
spherical  harmonics,  CfIr (9,  (p).  The  coordinates  9a  and  cpa  are  the  polar  and  azimuthal 
angles  corresponding  to  the  unpaired  B(2/;')  electron.  All  of  the  other  electronic  and 
nuclear  coordinates  have  been  integrated  over.  The  fonn  of  this  expansion  assumes  that 
mixing  of  states  outside  the  TP  manifold  will  be  weak.  Eq.  (  73  )  has  the  following  fonn 
when  represented  in  the  coupled  angular  momentum  basis: 
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(74) 
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The  inner  products  in  Eq.  (  74  )  are  evaluated  using  the  Wigner-Eckhart  Theorem 
yielding  Eqs.  (  75  )  and  (  76  ). 


C  (0,0) 
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=  (-l)"-s(2^  +  l)[(2;a  +  1)(2;„  +  1)]1/2  (  76  ) 
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W  s  f)\co  -n  ~(o')y 0  0  0/ 

In  Eqs.  (  75  )  and  (  76  )  3 -j  symbols  are  represented  as  (:::)  and  6 -j  symbols  are 

c  ^  37  82  83  91 

represented  as  ’  ’  ’  In  Appendix  E  the  symmetry  of  the  electrostatic  interaction 
potential  is  explored.  The  expansion  coefficients  Vxrxajl(r,  R)  are  found  by  expanding 
the  numerical  diabatic  PESs  Vxx,  Vzz,  Vyy,  and  Vxz  (  Eq.  (  37  )  )  in  terms  of  renonnalized 
spherical  hannonics  and  performing  a  tenn-by-term  comparison  with  the  analytic 
expansion  of  the  electrostatic  interaction  potential."  ’  ’ 

Comparing  the  full  Hamiltonian  as  expressed  in  Eq.  (  44  )  with  the  electronic 
Hamiltonian,  Eq.  (11),  the  following  connection  is  established: 


H, 


elec 


Te  +  Vel 

=  Hel2  +  V, 


+  fHz 

■  e 


+  Vel  +  B 


off 


el 


(77) 


As  explained  earlier,  Eq.  (  77  )  was  solved  using  ab  initio  techniques  resulting  in 
adiabatic  PESs.  The  ADT  mixed  these  surfaces  to  fonn  diabatic  PESs.  As  shown  in  Eq. 
(  37  ),  when  the  electronic  Hamiltonian  is  represented  in  the  diabatic  basis  is  has  the 
block  diagonal  form  shown  below. 
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In  order  to  make  a  comparison  to  the  numerical  PESs,  the  analytic  expansion  of  Eq.  (  77  ) 
in  the  asymptotic  basis  must  include  the  contributions  of  both  Vei  and  Hef.  This  is  done 
by  representing  the  electronic  Hamiltonian  using  the  asymptotic  electronic  basis  |S,  A)  = 

1 )  (x)  \{  —  1,  A  =  0,  ±1)  as  opposed  to  the  coupled  angular  momentum  basis  used  in 
Eqs.  (  74  ),  (  75  ),  and  (  76  ). 
v£A(r,e,R)  =  (S', A'\Vel  +  Hgi  |S, A) 

=  (S',  A' I  V  .|S,  A)  +  (S',A'|/?"2|S,A) 

j—1  \Qai  Q mj  | 


y  Vxrxai,(r,  R)C*r(9,  <p )  (A'| C%{Qa,  0„)|a)  8^  +  VH2(_r)8^8A>A 


(79) 


=  3X  X‘  i)(J  o“  l) 


+  kH2(r)5E'j;5A/A 

In  Eq.  (  79  )  the  results  of  Eq.  (  49  )  have  been  applied.  The  electronic  Hamiltonian  of 
the  hydrogen  molecule  becomes  the  ground  state  vibrational  potential  VH2  (r).  The 
electronic  part  of  the  electrostatic  interaction  potential  has  been  expanded  in  terms  of  3 -j 
symbols  using  the  Wigner-Eckhart  Theorem. 

The  potential  VH2  (r)  only  depends  on  the  hydrogen  bond  length  coordinate  r. 
Only  the  constant  renormalized  spherical  hannonic  term  Cq  ( 9 ,  0)  is  required  to  fit  this 
potential.  Eq.  (  79  )  can  be  written  as 
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The  8^it  terms  will  always  evaluate  to  unity  since  H2  is  restricted  to  its  ground  electronic 
state. 


The  properties  of  the  3 -j  symbols  lead  to  further  simplifications  to  Eq.  (  80  ).  The 
triangle  inequality  |1  —  Xa\  <  1  <  |1  +  Aa|  restricts  the  values  of  Aa  to  0,  1,  and  2.  3 -j 
symbols  also  have  the  property 


/ a  b  c\ 
VO  0  0/ 


0 


^0  0  0> 

when  the  sum  (a  +  b  +  c)  is  odd.  These  properties  lead  to 
l^A(r,0,7?)=^C^(0,0)  Varoo(r^)(-l)3_A'(3)1/2(_^  °  J)  Wa'a 


A' -A 


+  (r)Cg(#,  0)dA'A 


(81) 


(82) 


In  Eq.  (  82  ),  the  values  of  y  are  detennined  by  the  bottom  row  of  the  3 -j  symbol.  This 
property  requires  that  —A'  —  y  +  A  =  0,  or  /r  =  A  —  A',  where  A  and  A'  can  be  0  and  ±1 . 
Thus  y  takes  on  values  ranging  from  -2  to  2  in  integer  steps.  The  matrix  elements  of  the 
diabatic  electronic  potential  energy  A'  including  the  contribution  from  the  Tf?  molecule 
can  be  written  as 
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The  matrix  elements  of  VD  (R,  9,  r)  are  expressed  in  terms  of  coefficients  V\ali  given  by 

£<*•(0,0  =  0)^rf(r,«)  (84) 


Vka-/i 


To  cast  7°  (r,  0,  /?)  into  a  form  which  can  be  compared  with  the  numerical  ab  initio 
diabatic  PESs,  the  following  transformation  to  the  Cartesian  basis  is  made. 


|x>  =  —  (-id  +  i-D) 

|y>  =  -^(|l>  +  |-l» 

\z)  =  |0) 

This  transfonnation  yields  the  following  expression  for  VD  ( R ,  9,  r) : 
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As  expected,  VD  (r,  9,  R)  has  the  same  block  diagonal  form  as  Eq.  (  78  ).  Eq.  (  86  )  is  the 
analytical  fonn  of  the  electrostatic  Hamiltonian,  Eq.  (  77  ).  Now  the  numerical  diabatic 
PES  must  be  fit  to  the  same  functional  form  so  that  a  tenn-by-tenn  comparison  can  be 
made. 
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The  numerical  diabatic  PESs  are  expanded  in  terms  of  reduced  Wigner  rotation 
matrix  elements  which  have  the  following  relationship  to  renonnalized  spherical 
harmonics82, 83 

<5(R)  =  y4*(9'  <t>  =  °)  =  c4(e’  *  =  °)  (  87  ) 

The  following  functions  are  used  to  fit  each  diabatic  PES. 

Ezz(r,0,R)  =  £  l ^(r,R)d&(0) 

Ar=0 

Vxz(r,6,R)  =  7  k//(r,R)4(0) 

Ar=l 

14*  (r,  0,  R)  =  Es(r,  0,  R)  +  Pd(r,  0,  R) 

=  £  (E  R)  d& (0)  +  £  E/r  (r,  R)  (0)  (  88  } 

Ay = 0  — 2 


Eyy(r,  0,  R)  =  Vs(r,  9,  R)  -  Vd(r,  9,  R) 

=  ^  P/r(r,R)d^(0)-  £  v4(r,R)d^(9) 


The  functions  Ks. (r,  0,  R)  and  Ed(r,  0,  R)  are  introduced  to  separate  fitting  functions  with 
different  symmetry,  in  this  case  the  0^(0)  and  d^C 0)  terms.  They  are  defined  as 


1 40, 0,  R)  =  ^[Exx(E  9,  R)  +  Eyy(r,  0,  R)] 
Ed(r;  0,  R)  =  i[Pxx(r,  0,  R)  -  Pyy(r,  9,  R)] 
The  function  Vd(r,  0,  R)  can  also  be  defined  as 


(89) 


Vd(r,  9,  R)  =  \\VyyCr,  9,  R)  -  Vxx(r,  9,  R)]  (  90  ) 

Both  definitions  have  been  found  in  the  literature  (Weeks  uses  Eq.  (  89  )  while 

26  27 

Alexander"  ’  uses  Eq.  (  90  )  ),  and  the  choice  is  arbitrary. 
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The  fit  coefficients  in  Eq.  (  88  )  are  compared  term-by-term  with  the 
corresponding  expansion  coefficients  in  Eq.  (  86  )  yielding 

VzXzr(r,R)  =  VXr00(r,R)  +  ^VAr20(r,  R)  +  V„2(r)SAri0 

VXr(r,R)=^VAr21(r,R) 

1  (91) 
VXr(r,R)  =  VAr00(r,R)  - -VAr20(v,R)  +  E„2(r)5Ar,0 

Vdr<J-,R)=YV*r  22(T,r) 

These  equations  are  then  inverted  to  give  the  desired  fit  coefficients  VAr\ali(j,  R )  in  terms 
of  the  expansion  coefficients  of  the  numerical  PESs. 

VAr00(r,  R)=l  (vX/(r,  R )  +  2VXr(r,  R)')  -  VHi  (r)5A?.,0 

fAr20  (r,  R)=^  {vzz(j,  R )  -  V^Xr,  /?)) 

5  A  (92) 

VAr21(_r,R)=—v£(r,R) 

VAr2  2(r,R)=^=VdAr(r,R) 

As  seen  in  Eq.  (  92  )  the  definition  of  Vd(r,  0,  R )  introduces  a  sign  change  in  the 
fAr22  (a  R)  expansion  coefficients.  The  sign  of  VAr22  (r<  R)  was  verified  to  have  no  effect 
on  the  calculation  of  the  electrostatic  interaction  potential  Vel- 

Now  that  the  expansion  coefficients  VArxali  (r,  R)  are  known,  the  representation  of 
the  electrostatic  interaction  potential  Vei  in  the  coupled  angular  momentum  asymptotic 
basis  given  by  Eq.  (  74  )  is  complete. 
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Structure  of  the  Asymptotic  Representation 

The  full  Hamiltonian,  Eq.  (  60  ),  is  not  diagonal  when  expressed  as  a  matrix  in  the 
asymptotic  representation.  The  Kronecker  deltas  that  appear  in  the  matrix  elements  for 
each  term  in  the  full  Hamiltonian  detennine  which  asymptotic  basis  states,  defined  by  the 
labels  in  Eq.  (  57  ),  will  couple  with  one  another.  Terms  which  have  the  Kronecker  delta 
8X>X,  Eq.  (  63  ),  are  completely  diagonal  in  the  asymptotic  representation.  This  includes 
the  radial  and  vibrational  kinetic  energy  operators  (  Eqs.  (  62  )  and  (  64  )  ),  the  spin-orbit 
coupling  operator  (  Eq.  (  72  )  ),  the  rotational  kinetic  energy  operator  of  the  H2  molecule 
(  Eq.  (  65  )  ),  and  the  electronic  Hamiltonian  of  the  hydrogen  molecule  (  Eq.  (  79  )  )  are 
diagonal  in  the  asymptotic  representation. 

The  rotational  kinetic  energy  of  the  boron  atom  with  the  CM  of  the  H2  molecule 
(  Eq.  (71))  and  the  electrostatic  interaction  potential  (  Eq.  (  74  )  )  are  not  diagonal  in  the 
asymptotic  representation.  By  examining  the  Kronecker  deltas  in  these  terms,  the  matrix 
elements  can  be  organized  into  a  hierarchy  of  block  diagonal  pieces,  each  of  which  is 
infinitely  large.  The  /  block  are  the  largest  in  this  hierarchy.  Each  /  block  is  subdivided 
into  M  sub-blocks.  Each  M  block  contains  identical  infonnation  about  the  system,  so 
only  one  M  block  needs  to  be  considered.  This  is  illustrated  in  Figure  5  for  a  truncated  J 
=  +1/2  block.  The  grey  regions  in  Figure  5  are  zero.  Under  the  CS  approximation,  the 
gray  regions  within  a  given  M  block  are  zero;  however,  when  this  approximation  is  lifted 
the  entire  M  block  may  have  non-zero  values. 
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Figure  5.  A  visualization  of  the  7=  1/2  block  of  the  matrix  elements  of  the  asymptotic  Hamiltonian 

As  shown  in  Figure  5,  each  M  block  (red)  is  further  subdivided  into  smaller  P  blocks 
(green).  The  ellipses  indicate  that  each  P  block  is  infinite  in  dimension.  The  grey 
regions  indicate  the  location  of  non-zero  matrix  elements.  The  white  regions  are  where 
the  matrix  elements  are  zero.  The  electrostatic  interaction  potential  is  diagonal  in  this 
block.  The  Kronecker  delta  8P>P  in  Eq.  (  74  )  reinforces  a  result  that  arises  from  the  3 -j 
symbols  in  Eq.  (  75  )  and  (  76  ).  The  following  two  equations  are  obtained  using  the 
condition  that  the  bottom  row  of  a  3 -j  symbol  must  sum  to  zero: 

—k'  +  fi  +  k  =  0  ,  , 

a)-n-o)'  =  0  (  ’ 

These  equations  imply  that  k'  +  co'  —  k  +  co,  or  A P  —  0.91  However,  the  matrix 
elements  of  the  tumbling  motion  of  boron  and  H2  molecule  are  not  diagonal  in  the  P 
block.  This  is  the  primary  reason  for  invoking  the  CS  approximation.  Under  the  CS 
approximation  off-diagonal  P  blocks  are  zero,  and  a  specific  P  block  can  be  chosen  for 
the  calculation  rather  than  having  to  consider  the  entire  M  block. 
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The  P  block  can  be  further  subdivided  into  a  parahydrogen  block  and  an 
orthohydrogen  block.  Since  these  states  do  not  couple  they  are  block  diagonal  in  P.  In 
this  work  only  transitions  from  the  ground  rotational  state  j  =  0  re  considered.  This 
restricts  the  basis  functions  to  include  only  parahydrogen  states  (J  =  0,2,4, 6  ...). 

Given  that  only  even  values  of j  will  be  used  and  that  ja  will  be  either  1/2  or  3/2, 
the  basis  functions  for  the  J  =  1/2  block  under  the  CS  approximation  can  be  identified. 
This  work  examines  the  P  =  +1/2  block,  which  leads  to  the  requirement  that  k  +  a>  = 

+  1/2.  The  basis  functions  which  meet  this  requirement  are  shown  below. 


7=0, 
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As  shown  in  Figure  6,  the  parahydrogen  block  has  a  2-by-2  sub-block  corresponding 
to  j  —  0.  For  each  additional  value  of  j  considered  an  additional  six  basis  functions  are 
added.  At  first  glance,  it  does  not  seem  that  some  of  the  basis  functions  in  Eq.  (  94  )  and 
(  95  )  belong  in  the  J=  1/2  block  since  their  values  for  j  and  ja  sum  to  greater  than  1/2. 
Their  inclusion  in  this  block  is  a  consequence  of  how  the  three  angular  momenta  of  this 

83 

system  couple  to  give  the  total  angular  momentum. 
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Figure  6.  A  visualization  of  the  parahydrogen  block  within  the  P  sub-block  of  the  matrix  elements  of  the 
asymptotic  Hamiltonian 


Since  an  infinite  number  of  basis  function  cannot  be  considered  in  a  numerical 
calculation,  the  basis  set  must  be  truncated.  This  is  done  by  considering  values  of  j  for 
which  the  rotational  energy  is  less  than  the  total  energy  of  the  collision  being  considered. 
This  threshold  value  of  j  is  designated  jmax  and  leads  to  a  basis  size  given  by 
d  —  2  +  3  jmax.  The  rotational  energy  levels  (in  atomic  units)  of  the  H2  rotor  are  given 


by 


E.(r)=j(j  +  1) 

]{  J  2 Mh  r2 


(96) 


These  energies  depend  on  the  bond  length  of  Tb  which  is  no  longer  a  constant.  However, 
for  the  energies  considered  in  this  work,  the  equilibrium  bond  length  is  used  as  an 
estimate.  Energies  between  0.0  a.u.  and  0.01  a.u.  are  considered  (less  than  the  energy 
required  to  excite  the  first  vibrational  mode  of  H2).  Using  Eq.  (  96  ),  this  leads  to 
Jmax  =  6  and  a  basis  size  of  20. 
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The  Time-Dependent  Schrodinger  Equation 


The  time-dependent  Schrodinger  equation  (TDSE)  can  be  represented  using  the 
basis  functions  defined  by  Eqs.  (  94  )  and  (  95  ).  The  general  form  of  the  TDSE  is  given 
by  Eq.  ( 2  ).  The  Hamiltonian  can  be  broken  into  a  kinetic  energy  tenn  and  a  potential 
energy  term.  For  the  B  +  H2  system  the  kinetic  term  is  composed  of  the  kinetic  energy 
operators  given  by  Eqs.  (  62  )  and  (  64  ).  The  potential  energy  term  is  the  sum  of  the  H2 
rotational  and  tumbling  kinetic  energy,  the  electrostatic  interaction  potential,  and  the 
spin-orbit  coupling.  The  matrix  elements  of  this  operator  are  referred  to  as  effective  PES 
and  are  given  by 


VeDff(r,R)a  =  U'.r'.R' 


L2  f 

+  - - -  +  Vei  +  Hgf  +  Hs 


C,  r,R 


(97) 


2  [ib,h2R2  2[iH2r2 

The  labels  f  and  f '  represent  states  for  which  appropriate  values  for  each  of  the  labels  in 
Eq.  (  57  )  have  been  chosen.  The  superscript  D  is  used  to  indicate  that 
not  diagonal  in  this  representation. 

In  matrix  form,  the  TDSE  is  given  by  (atomic  units) 


a  AtqC r,R,t)\  ///„  H12 

®2(EP,t) J  =  (  ^2i  R22  - j(  <P2(r,R,t)J 

(  -1  d2  -1  32}/1  0 

(2  ^B,H2dR2  2^H2dr2\y°  1  .) 

/Pe//(r<P)ll  Pe//(r)P)l2 
+  \Veff(r>R)  21  Vgff(r,R)  22 


\]  R,  t)\ 
I  ( <t>2(r,R,t)  J 


(98) 


where  the  radial  kinetic  energy  and  the  vibrational  kinetic  energy  operators  have  been 
simplified  by  introducing  the  reduced  wave  function  0(r,  R,  t )  =  rJR'P(r,  R,  t).  The 
formal  solution  (in  atomic  units)  governing  the  time  evolution  of  the  reduced  wave 
functions  is 
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(99) 
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Eq.  (  99  )  represents  the  time  evolution  of  a  collection  of  wave  functions,  denoted 
Oj(r,  R,  At ),  under  the  influence  of  the  Hamiltonian  H.  The  time  interval  is  given  by  At. 

The  wave  functions  propagate  on  a  potential  surface  given  by  a  corresponding 
diagonal  element  of  V(?^-  (r,  R).  However,  with  each  propagation  step  At  the  off- 
diagonal  elements  of  Vg^(r,  R)  couple  the  wave  functions  causing  probability  amplitude 
to  be  redistributed  over  the  entire  collection  of  PESs.  If  the  effective  potential  energy 
Vgff  (r,  R)  were  diagonal,  this  redistribution  would  not  occur  and  each  wave  function 
would  propagate  independently  on  its  own  potential  surface. 

Eq.  (  99  )  can  be  expressed  in  a  form  more  suitable  for  numeric  computation  by 
applying  the  split  operator  approximation  (SOA). 
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oAr.fl,  o)\  (100) 
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In  Eq.  (  100  ),  tenns  of  order  than  (At)3  and  higher  involve  commutators  off  and  V. 
These  terms  can  be  neglected  provided  that  the  time  step  At  is  sufficiently  small.  The 
SOA  is  described  in  detail  in  Leforestier  et  al.  Wave  packet  propagation  using  the 
SOA  has  been  demonstrated  by  Alvarellos  and  Metiu. 

Each  exponential  function  in  Eq.  (  100  )  is  easily  evaluated  by  moving  into  the 
representation  where  each  operator  is  diagonal.  An  exponential  function  with  a  matrix 
argument  must  be  expanded  in  a  Taylor  series  to  evaluate  it.  If  the  matrix  is  diagonal  this 
Taylor  series  simply  yields  a  new  diagonal  matrix  where  the  values  of  each  element  are 
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the  exponentiated  values  of  the  corresponding  value  of  the  original  matrix.  This  is 
illustrated  in  Eq.  (  101  ). 
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The  kinetic  energy  operators  in  Eq.  (  98  )  are  diagonal  in  the  asymptotic  representation. 
A  Fourier  transfonn  changes  the  kinetic  energy  operators  for  the  coordinate 
representation  into  the  momentum  representation,  replacing  the  derivative  operators  — V2 
with  numbers  k2.  This  simplifies  the  evaluation  of  the  exponential  function  by 
multiplying  momenta  rather  than  differentiating  with  respect  to  nuclear  coordinates. 

The  exponential  functions  containing  the  effective  potential  energy  term  can  also 
be  transformed  in  a  similar  way.  It  is  possible  to  define  a  unitary  transformation  U  for 
each  (r,  R )  point  that  diagonalizes  the  effective  potential  energy.  Using  this  spatially 
dependent  transformation  matrix  U(r,  R),  the  effective  potential  Ve//,  Eq.  (  97  ), 
transform  as  follows 

U+(r,ff)VeD//<I>D  =  U+(r,R)VeD//U(r,R)U+(r,R)<I>D 


(102) 


Once  the  effective  potential  has  been  transfonned,  Eq.  (  101 )  can  be  used  to  evaluate  the 
exponential. 

The  time  evolution  of  the  reduced  wave  functions  under  the  SOA  is  given  by 
U»D(r,  R,  At)  =  UUte“^At/2UlU+7:'_17:'e_ifAt7:'“1T'UlUte“^A£/2lUU+<l>D(r,  R,  0)  (  103  ) 


Eq.  (  103  )  describes  each  processing  step  as  4>°(r,  R,  0)  becomes  4>°(r,  R,  At).  First, 
the  wave  packet  and  exponential  operator  e~lVAt/2  are  transfonned  to  the  adiabatic  basis 
and  then  multiplied.  The  result  of  this  operation  is  first  transformed  back  to  the  diabatic 
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representation  and  then  Fourier  transformed  to  the  momentum  representation.  A  two 
dimensional  Fourier  transform  used  to  transform  both  radial  and  vibrational  kinetic 
energies,  Eqs.  (  62  )  and  (  64  ).  The  exponential  operator  e~lTAt  is  diagonalized  by  this 
transformation.  An  inverse  Fourier  transfonn  returns  the  result  of  this  multiplication 
back  to  the  coordinate  representation.  Here  again  the  wave  packet  and  exponential 
operator  e~lVAt /2  are  transfonned  to  the  adiabatic  basis  and  then  multiplied.  Finally,  the 
result  is  transformed  back  into  the  diabatic  representation  yielding  the  time  evolved  wave 
packet  <l>D(r,  R,  At) . 

Introduction  to  the  Scattering  Problem 

By  using  the  propagation  scheme  described  by  Eq.  (  100  ),  an  initial  wave  packet 
can  be  propagated  under  the  influence  of  the  effective  potential,  Eq.  (  97  ).  In  this  way 
the  nonadiabatic  effects  captured  in  the  full  Hamiltonian  are  incorporated  into  the  wave 
packet  dynamics.  Infonnation  about  the  system  can  be  extracted  from  these  evolved 
wave  packets.  This  work  focuses  on  calculating  scattering  matrix  (S-matrix)  elements. 
S-matrix  elements  are  used  to  calculate  the  probability  of  state  transitions  due  to  an 
interaction.  The  method  of  calculating  S-matrix  elements  from  the  propagation  of  wave 
packets  falls  within  the  framework  of  the  CPM.  Before  discussing  the  details  of  the 
CPM,  a  very  brief  introduction  to  scattering  theory  is  given.  All  of  the  scattering  theory 
presented  here  is  discussed  in  more  detail  in  Taylor.94 

In  atomic/molecular  scattering  experiments  the  interaction  between  reactants 
takes  place  within  a  region  no  larger  than  several  atomic  diameters  and  on  very  short 
time-scales.  Although  new  ultra-fast  spectroscopic  techniques  are  beginning  to  probe 
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these  small  time-scales,  nearly  all  of  what  happens  within  the  interaction  region  is 
unobservable.  The  fundamental  goal  of  a  scattering  problem  is  to  relate  the  state  of  a 
system’s  reactants  well  before  the  scattering  event  to  the  state  of  a  system’s  products  well 
after  the  interaction.  By  doing  so,  the  initial  state  can  be  related  to  the  possible 
experimental  outcomes  when  the  final  state  is  measured. 

Figure  7  depicts  a  wave  function  as  it  propagates  through  the  interaction  region 
(shown  as  a  grey  box). 

A- 

Figure  7.  A  visualization  of  a  wave  packet  as  it  scatters  off  of  the  interaction  region 

The  time  t  ->  —  oo  refers  to  a  time  before  the  interaction  when  the  wave  packet  is  spatially 
separated  from  the  interaction  region  (has  no  significant  probability  amplitude  in  the 
region).  The  time  t  —>  +co  refers  to  a  time  after  the  interaction  when  the  wave  packet  has 
left  the  interaction  region.  Both  of  these  regions  are  referred  to  as  the  asymptotic  limit. 

In  the  asymptotic  limit  the  wave  packet  evolves  under  the  influence  of  the  asymptotic 
Hamiltonian  denoted  by  H0 .  The  asymptotic  Hamiltonian  is  the  limit  of  the  full 
Hamiltonian  when  the  wave  packet  is  spatially  separated  from  the  interaction.  By  the 
time  the  wave  packet  enters  the  shaded  interaction  region  of  Figure  7  the  full  Hamiltonian 
H  has  lost  its  asymptotic  character.  The  choice  of  the  time  origin  t  —  0  is  arbitrary.  A 
common  choice  places  the  wave  packet  well  within  the  interaction  region  at  t  =  0.  The 


|<x|'tJ(t  =  0))|2  |(x|4J(t  ->  +oo))|2 
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wave  packet  continues  to  evolve  as  it  exits  the  interaction  region.  As  the  wave  packet  is 
spatially  separated  from  the  interaction  region  the  full  Hamiltonian  again  regains  its 
asymptotic  character.  At  some  time  well  after  the  wave  packet  has  exited  the  interaction 
region  (t  -»  +oo)  it  is  compared  with  a  product  state  wave  packet. 

Figure  7  depicts  a  single  channel  scattering  event — a  single  input  wave  packet 
scattering  into  a  single  output  packet.  The  B  +  H2  system  has  many  different  channels. 
The  number  of  channels  available  is  determined  by  the  value  of  jmax.  For  instance,  if 


imax  —  0  then  there  are  two  possible  channels  to  consider: 


.  Figure  8 


depicts  a  two  channel  scattering  event  where  the  incoming  wave  packet  is  defined  in  a 
single  channel. 


|<x|'F2(t-»  +°o))| 
\ 


Kxft'At  ->  +oo)>| 


Figure  8.  A  visualization  of  a  two  channel  scattering  event 


2 


2 


As  the  wave  packet  enters  the  interaction  region,  the  probability  amplitude  of  the  wave 
packet  is  redistributed  among  both  channels  as  time  progresses.  Both  wave  packets  will 
continue  to  propagate  giving  a  detector  a  chance  of  measuring  an  outcome  for  either 
channel.  For  jmax  =  6  there  are  20  different  channels  (when  no  vibrational  transition 
occurs,  v  =  v')  each  labeled  by  a  different  set  of  quantum  numbers  contained  in  f 
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(  Eq.  (  57  )  ).  The  wave  packets  shown  in  Figure  7  and  Figure  8  are  one-dimensional. 

The  B  +  FE  system  has  two  dimensional  wave  packets.  The  dimensions  correspond  to 
the  Jacobi  coordinates  R  and  r. 

The  reactant  states  of  the  B  +  H2  system  are  defined  by  forming  a  superposition  of 
eigenstates  of  the  asymptotic  Hamiltonian  (  Eq.  (  58  )  ).  The  reactant  state  can  be 
specified;  however,  the  product  state  is  a  result  of  the  scattering  event.  The  reactant  and 
product  states  are  defined  as 


l^react)  =  |V(t  ->  -oo)> 

I'kprod)  =  \^(X  "»  +°°)> 

The  TDSE,  Eq.  (  98  ),  prescribes  how  to  propagate  this  wave  packet  in  time.  In  Eq. 


(104) 


(  99  ),  the  operator  e~lHAt  is  called  the  time  evolution  operator  (TEO).  As  shown  in  Eq. 

(  99  )  the  TEO  propagates  a  wave  function  defined  at  t  =  0  to  a  time  t  =  At.  The  TEO  is 
used  to  propagate  the  reactant  state  through  the  interaction  region.  The  reactant  state  will 
eventually  evolve  into  the  product  state. 

The  reactant  wave  packet  is  defined  as  a  superposition  of  plane  waves  labeled  by 
k  within  a  given  channel  defined  by  and  v  given  by  Eq.  (  105  ). 


react 


=  J  d/c^(/c^)|£,  v^,  «T) 


<E> 


react 


(105) 


The  expansion  coefficients  of  this  superposition  of  plane  waves  are  given  by  t]y{k^  = 
(^^Keact)- 
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The  discussion  of  scattering  in  Taylor94  defines  an  initial  wave  packet  |cpin)  at 
t  =  0  that  when  propagated  to  the  asymptotic  limit  using  the  asymptotic  Hamiltonian 
produces  the  reactant  state  of  system  |V(t  —  oo)).  In  the  asymptotic  region 
| (1)^ (t  ~ °°))  —  |^(t  ->  — oo)).  Choosing  an  initial  wave  packet  already  localized  in 
the  asymptotic  region  eliminates  this  step.  The  TEO  is  then  used  to  propagate 
|  chin  (t  ~ °°))  forward  in  time  under  the  influence  of  the  full  Hamiltonian  to  create  the 
actual  state  of  the  system  during  the  interaction  ^(t  =  0)).  This  process  is  depicted  in 
Figure  9. 


Figure  9.  A  depiction  of  how  the  actual  state  of  the  system  |cp(t  =  0))  is  mapped  to  |<t>in) 

When  grouped  together,  the  TEOs  that  accomplish  the  propagations  depicted  in 
Figure  9  are  called  Moller  operators.  The  TEO  are  given  separately  by  Eqs.  (  106  )  and 
(107  ). 

U0(t)  =  e~iRot  (106) 

Eq.  (  106  )  represents  the  TEO  under  the  asymptotic  Hamiltonian.  The  TEO  under  the 

full  Hamiltonian  is  given  by  Eq.  (  107  ). 

U(t)  =  e~iRt  (  107  ) 

Using  Eqs.  (  106  )  and  (  107  )  the  Moller  operators  are  defined  as  follows. 
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(108) 


fl+ =  lim  U^(t)U0(t) 

-  £->+co 

The  M0ller  operator  that  maps  |  O  j  n )  to  |T(t  =  0))  is  H+.  The  sign  in  the  subscript  of  the 
Moller  operator  is  opposite  the  sign  in  the  limit.  The  Hennitian  conjugate  of  the  full 
Hamiltonian  TEO  is  used  to  ensure  that  it  propagates  in  the  opposite  time  direction  as  the 
asymptotic  Hamiltonian  TEO.  The  state  resulting  from  the  operation  of  fl+  on  |Oin)  is 
called  the  Moller  state  |4>+). 


|0+)  =  n+|cpin>  (109) 

The  Moller  state  |d>+)  represents  |V(t  =  0)),  the  actual  state  of  the  system  at  t  =  0. 

Now  the  Moller  state  must  be  propagated  forward  in  time  to  the  asymptotic  limit 
to  determine  the  final  state  of  the  system  |V(t  ->  +oo))  which  is  the  product  state 
l^prod)-  This  final  state  ^(t  -»  +oo))  must  also  be  mapped  to  a  wave  function  defined 
at  t  =  0.  Both  of  these  steps  are  done  by  applying  the  Moller  operator  fit  on  the  state 
|<f>+)  to  give  |0out).  The  state  |<bout)  is  the  state  defined  at  t  =  0  that  when  propagated 
forward  in  time  under  the  asymptotic  Hamiltonian  becomes  the  product  state. 

Eq.  (110)  shows  how  the  Hennitian  conjugate  of  the  Moller  operator  f l_ 
switches  the  order  of  propagation. 


at  =  j  lim  f/T(t)  f/0(t)j 
=  lim  [/J(t)t/(t) 

t->  +  °0 


(110) 


Now  the  wave  packet  is  propagated  forward  in  time  under  the  full  Hamiltonian  and  then 
propagated  backward  in  time  using  the  asymptotic  Hamiltonian.  Thus  the  Moller 
operator  Hi  maps  |<f>+)  onto  |<t>out).  Figure  10  depicts  how  this  happens  for  |<3>out).  The 
states  |Oin  oUt)  represent  wave  functions  that  when  propagated  to  their  respective 
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asymptotic  limits  map  onto  the  asymptotic  wave  function  of  the  actual  state  of  the  system 
^(t  -*  +oo)).  These  states  are  often  referred  to  as  the  ‘in’  or  ‘out’  asymptote. 


Figure  10.  A  depiction  of  how  the  actual  state  of  the  system  |'F(t  =  0))  is  mapped  to  |<t>out) 

The  scattering  operator  S  is  constructed  by  combining  both  Moller  operators. 

|chout)  =  ntn+|oin) 

=  $  l*in>  (  j 

The  scattering  operator  S  maps  the  ‘in’  asymptote  |Oin)  to  the  ‘out’  asymptote  |$>out). 

Eq.  (Ill)  can  be  used  to  calculate  the  probability  of  measuring  an  outcome  of  a 
scattering  event.  This  probability  is  expressed  as 

p  =  iKishn>i2 

=  l<<bml<bout>l2  (112) 

=  |((h_|0+)|2 

In  Eq.  (112),  the  state  |d>m)  is  the  experimental  outcome  state.  It  is  chosen  from  the  set 
of  eiegnfunctions  of  a  Hennitian  operator  (an  observable).  For  example,  if  an  experiment 
measured  the  total  energy  of  the  final  state  after  the  interaction,  the  state  |Om)  would  be 
chosen  from  the  eigenfunctions  of  the  Hamiltonian  operator.  The  initial  state  |<3>in)  and 
the  measurement  state  |Om)  can  also  be  chosen  as  an  eiegenstate  of  the  asymptotic 
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Hamiltonian.  This  would  give  the  probability  of  system  starting  in  the  state  being 
measured  in  the  state  |<J>m)  after  the  interaction. 


The  Time-Dependent  Channel  Packet  Method 

In  the  time-dependent  CPM,  reactant  states  (the  channel  packets)  are  chosen  and 
propagated  using  TEOs. 43,44  These  evolved  wave  packets  are  then  compared  to 
measurement  wave  packets  defined  in  many  channels  to  calculate  probabilities  using  Eq. 
(  112).  The  basis  functions  given  by  Eqs.  (  94  )  and  (  95  )  define  the  input  and  output 
channels  of  the  scattering  problem  for  the  B  +  H2  scattering  problem.  Wave  packets  are 
created  within  these  channels. 
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A  superposition  of  plane  waves  is  required  to  form  a  normalizable  wave  function.  In 
conjunction  with  the  reduced  wave  function  used  in  Eq.  (  98  ),  the  definition  of  the 
reduced  momentum  operator  pR  in  both  the  momentum  and  coordinate  representations  is 
given  as  (in  atomic  units) 


Pr  kR  (momentum) 
pR  <->  —  i^-(  coordinate) 

uR. 


(115) 


The  action  of  the  asymptotic  Hamiltonian  can  then  be  expressed  as 

#o|<>C,kC)  =  (Hq  +  Hint)\(,V^,K^) 

=  (4  +  4t)  If.vV) 


(116) 


where  E^nt  is  the  energy  eigenvalue  associated  with  Hint.  The  energy  eigenvalue 
associated  with  Hq  is  £j< .  The  label  (  is  retained  on  the  eigenvalue  because  it  is 
associated  with  a  specific  channel  given  by  |<(,  ,  k 

When  vibrational  transitions  are  permitted,  the  label  creates  manifolds  within 
the  channel  thus  multiplying  the  size  of  the  basis  by  the  number  of  vibrational 
transitions  considered.  This  work  will  restrict  H2  to  the  ground  vibrational  state  (v^  = 
v *>'  —  0).  However,  the  theory  developed  thus  far  is  general  enough  to  accommodate  the 
calculation  of  transitions  to  higher  order  vibrational  modes.  In  these  cases  the  label 
designates  a  channel  and  designates  a  manifold  within  the  channel.  There  is  a  risk  of 
entering  the  energy  regime  where  boron  and  H2  react  to  fonn  BH2  when  modeling  higher 
order  vibrational  transitions.  Modeling  the  reaction  requires  additional  potential  energy 
surfaces  not  considered  in  this  work. 

Each  plane  wave  component  is  associated  with  an  energy  given  in  atomic  units  by 
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(117) 


E  =  E(  |  l*Sl 
int  2EBjH2 


,1/2 


\^\={2^H2(e-eQ} 

Fonning  a  superposition  of  plane  waves  serves  two  purposes:  it  creates  a  normalizable 
wave  packet,  and  allows  the  calculation  of  scattering  probabilities  as  a  function  of 

energy.  The  initial  and  outcome  wave  packets  |o[n  mj  are  constructed  by  fonning  a 

linear  combination  of  plane  waves.  Using  eigenfuntions  defined  at  t  =  0  this 
superposition  can  be  expressed  as 


where 


)  =  j  dK^tV^K^)l^y,K 
=  j  dK^l(K^,V^,K^) 

=  j  dE&±  (4)|c,^,4 


clT 


in,m 


(118) 


(119) 


and 


=  krf.K 


<t>; 


111,111 


ril  (4)  =  t,v^,4 


in,m 


(120) 


’’iW = (!w)2’’hK0 


Eq.  (118)  expresses  the  superposition  as  a  function  of  the  wave  number  k ^  and  the 
energy  E Eqs.  (119)  and  (  120  )  give  the  relationships  needed  to  perfonn  this 
transformation.  Here  the  label  (  has  been  retained  in  the  superscript  for  |o[nm) 
indicating  that  the  wave  packet  is  defined  in  a  specific  channel.  The  expansion 
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coefficients  in  Eq.  (  120  )  are  chosen  so  that  77+ (k1^)  corresponds  to  |d>[n)  and  17^  (k^) 
corresponds  to  The  coefficients  77+ (jc1^)  are  generally  chosen  to  give  a  Gaussian 

function  in  the  coordinate  representation. 

Once  these  states  have  been  constructed,  they  can  be  propagated  using  the  Moller 
operators.  Theoretically  the  propagation  must  be  done  over  an  infinite  time  interval; 
however,  for  computational  purposes  the  wave  packets  must  be  created  in  and  propagated 
well  into  the  regime  of  the  asymptotic  Hamiltonian.  For  the  B  +  H2  system,  this  occurs 
as  the  electrostatic  interaction  potential  approaches  zero.  Once  propagated,  Eq.  (112) 
can  be  used  to  calculate  the  probability.  However,  this  is  not  the  probability  of  a  single 

eigenstate  scattering  into  a  different  eigenstate.  The  states  |(f»[nm|  are  formed  from  a 

superposition  of  eigenstates  of  the  asymptotic  Hamiltonian.  It  is  possible  to  extract 
energy  (momentum)  resolved  S-matrix  elements — the  probability  of  a  well-defined 
energy  state  transitioning  to  another  well-defined  energy  state. 

Energy  Resolved  S-matrix  Elements 

The  eigenstates  of  the  translation  operator  (  Eq.  (  114))  are  not  proper  wave 
functions  since  they  cannot  be  nonnalized.  As  a  result,  a  superposition  of  these  states 
was  formed  to  create  a  normalizable  wave  packet.  The  wave  packet  itself  is  not  an 
eigenstate  of  the  asymptotic  Hamiltonian  and  does  not  have  a  well  defined  energy.  In 
this  section  the  method  for  extracting  energy  (momentum)  resolved  S-matrix  elements  is 
derived.  The  discussion  in  this  section  is  motivated  by  the  derivation  of  energy  resolved 
S-matrix  elements  performed  by  Tannor  and  Weeks  1993. 42 
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The  derivation  begins  by  examining  the  different  propagation  steps  contained  in 


the  scattering  matrix  operator  in  the  term  S  o[nj  found  in  Eq.  (112).  The 

superscripts  f  and  f '  serve  as  a  reminder  that  the  initial  state  and  the  measurement  state 
can  be  defined  in  different  reactant/product  channels.  These  superscripts  will  be  omitted 
to  simplify  the  notation.  To  simplify  the  notation  further  |Oa)  will  represent  |o[n)  and 

|  Op)  will  represent  |  O^V 


(Op|5|oa)  =  (Op|ntn+|oa) 

=  (Op|fltfi+|Oa) 

=  H  (,!!"„  PjCOPftOjflim  UHt)  Soft)}  |<fa) 

=  (opCt  -  +»)|  film,,  Pa')}  film  pt(t)}  Ka  -> 

=  /op(t  -*  +oo)  |  lim  Oa(t  — >  —  oo)| 


(121) 


Eq.  (  121  )  shows  how  the  different  propagation  steps  of  the  Moller  operators  can  be 
separated  and  used  to  propagate  the  wave  packets  |Oa)  and  |Op)  to  their  respective 
asymptotic  limits.  The  remaining  TEOs  can  be  consolidated.  The  limits  on  t  and  t' 
create  a  positive  time  interval,  the  propagation  time  tp .  The  propagation  time  is  given  by 

tv  =  t'-t  (  122  ) 

With  this  new  definition  of  time,  Eq.  (  121  )  can  be  written  as 

=  t  hmJOp(t  ->  +oo)|e-lHtp|cpa(t  ->  -oo)) 

"  r  (123) 

where  the  time-dependent  correlation  function  C^Atp)  is  defined  as 
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cc'c(tp)  =  (*p(f  +°°)|e  lHtp\<S>a(t  -»  -00))  (124) 

The  correlation  function  C^?(tp)  captures  the  overlap  between  the  measurement  state 

|$p(t  +  oo))  and  the  evolving  reactant  state  |(h^eact}  =  | (t  ->  — oo))  as  a  function 

of  time.  When  the  limit  as  the  propagation  time  tp  — >  +oo  is  taken,  the  result  is  the  S- 

matrix  element  (Op|s|Oa). 

The  derivation  proceeds  by  taking  the  Fourier  transform  of  both  sides  of  Eq. 

(  121  )  arriving  at  the  following  expression 


lim 

tp->  +  0°  J 


dt'v{  <P^\S\<Pa)eiEtv 


lim 

„->  +  00 


j 


dtp  C, 


(125) 


o  o 

The  time  origin  is  chosen  so  that  the  propagation  of  |Oa(t  ->  —  oo))  begins  at  tp  —  0. 
However,  the  limits  of  integration  can  be  extended  to  earlier  times  since  there  is  no 
significant  overlap  between  |<t»p(t  ->  +oo))  and  the  evolving  reactant  state 
|Oc((t  ->  —  oo))  for  times  earlier  than  tp  —  0.  Thus,  Eq.  (  125  )  can  be  cast  into  a  fonn 
easily  recognized  as  a  Fourier  transform. 

fp  ip 


lim 

tp->  +  oo 


dt!p{'Pfi\S\'Pa)eiEtv  = 


lim 

tp->+  00 


(126) 


dtp(<P^\S\^a)eiEtv  =  dtpC^^(tp)  eiEtv 


The  right-hand  side  of  Eq.  (  126  )  is  the  Fourier  transform  of  the  correlation  function. 
This  can  be  computed  numerically  by  propagating  the  reactant  wave  packet  using  the 
TEO  e~lHtP,  calculating  the  correlation  function  C^>^(tp)  for  each  value  of  tp,  and  then 


62 


taking  the  finite  Fourier  transfonn  of  the  resulting  array  of  points.  The  right-hand  side  of 
Eq.  (  126  )  can  be  manipulated  further  to  reveal  energy-resolved  S-matrix  elements. 

This  is  done  by  first  expressing  (d>p|s|d>a)  in  terms  of  the  asymptotic  eigenstates 
by  using  completeness  relationships. 


OO  00 


(<&p|$|<Da)=  £  j  dsfj 


z'.P'v*  ,vr  0 

x  (a',va'  ,Eg' |Oa) 
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The  tenn  (/?',  v^' ,  S  a',  va  ,  Eg  )  is  the  energy-resolved  S-matrix  element.  The  terms 


(a',  va' ,  Eg'  |d>a)  in  Eq.  (  127  )  can  be  express  as  follow  by  using  the  definition  given  by 
Eq.  (  118). 
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In  a  similar  fashion,  the  tenn  /<f>p  (t  — >  +oo)  /?',  ,E^  \  can  be  written 


(*4*  -  +oo)|/?',v*>f )  =  r,p+  (E?)  8p'p8vP'vft  (  129  ) 

When  Eqs.  (  128  )  and  (  129  )  are  inserted  into  Eq.  (  127  )  the  following  expression  is 
obtained. 
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Eq.  (  130  )  is  inserted  in  the  right-hand  expression  of  Eq.  (  126  )  to  yield 
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(131) 
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The  term  E  can  be  defined  in  the  following  way 


E  =  E"  -  El 


(132) 

where  E"  is  a  constant  energy  offset.  The  integral  over  time  now  leads  to  a  delta  function 
as  shown  in  Eq.  (  133  ). 
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The  delta  function  S{E"  —  E£  )  leads  to  the  following  simplification. 
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Eq.  (  135  )  introduce  the  on-shell  S-matrix  element  Sap(E)  defined  as 


Sap{E)s(E^  -E)  =  ^',vP',EZ'\s\a',va',E)  (  135  ) 

The  on-shell  S-matrix  element  allows  for  further  simplification  of  Eq.  (  134  ).  The  right- 
hand  term  of  Eq.  (  126  )  becomes 


I  dt^\S\<Pa)eiEtv  =  2tci)Z (E)V“ (E)Sap(E) 

—  00 

Further  manipulation  of  Eq.  (  126  )  yields 

+  00 

J  dtpC^^(tp)  elEtP 


(136) 


■Sa/?(E)  — 


2nriZ(E)ri“(E) 


(137) 


Eq.  (  137  )  can  also  be  expressed  in  terms  of  momentum  by  using  Eq.  (  120  ). 
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sections  are  given  by 

=  V^yYj^21  +  l5«'(£,)|2  (  139  ) 

; 

where  is  defined  in  Eq.  (  117).  The  total  angular  momentum  J  is  contained  within 
the  state  labels  (  and  ("  defining  the  asymptotic  basis  set.  Reaction  rates  can  be 
calculated  using  the  following  equation. 

00 

/  8 k  T  \  C  ~E 

k^{T)  =  (  -  J  (kBT)~2  EeWa  ,{E)dE  (  140) 

\nE-B,H2J  J 

0 

This  work  restricts  the  total  angular  momentum  to  J  —  1/2.  Consequently  cross-sections 
and  reaction  rates  will  not  be  calculated. 
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III.  Results  and  Discussion 


The  Adiabatic  Potential  Energy  Surface  and  Derivative  Coupling  Data  Set 

The  1  A 2  . A and  1  ~A  "  adiabatic  PESs  and  DCTs  used  in  this  research  were 
calculated  by  Dr.  Dave  Yarkony  (Johns  Hopkins  University)  using  his  suite  of  ab  initio 
programs.  The  adiabatic  PESs  were  calculated  at  the  SA-MCSCF/CI  level.  The 
contracted  basis  set  (13s8/?3t/)/[8s5/?3<i]  was  used  for  boron.95’96  The  contracted  basis  set 
(%s2>p\d)/[6s2>p\d\  was  used  for  hydrogen.  This  approach  for  the  electronic  structure 
calculation  was  used  in  previous  work  on  the  B  +  H2  system  and  has  been  validated  by 
experimental  endoergicity  measurements.29’ 45  Previous  electronic  structure  calculations 
were  carried  out  on  a  two  dimensional  grid  given  by  Jacobi  coordinates  ( R ,  9);  the 
coordinate  r  was  fixed  at  the  equilibrium  value  of  the  H2  bond  length.  ’  ’  The 
electronic  structure  calculation  for  this  work  was  carried  out  for  several  values  of  the  H2 
bond  length  making  it  possible  to  consider  H2  vibrational  transitions  in  future  work. 

Adiabatic  PESs  and  DCTs  were  calculated  at  points  on  a  three  dimensional  grid 
given  by  Jacobi  coordinates  (r,  6,  R ).  Each  point  on  this  grid  represents  a  given  spatial 
configuration  of  the  B  +  H2  nuclei.  The  values  of  R  range  from  zero  to  10  a.u.  in  steps  of 
0.2  a.u.  for  a  total  of  5 1  grid  points  in  this  coordinate.  The  values  of  9  range  from  zero  to 
n/2  radians  in  steps  of  n/20  for  a  total  of  1 1  grid  points  in  this  coordinate.  The  values 
of  r  include  0.7  a.u.  and  range  from  0.904  a.u.  to  4.402  a.u.  in  steps  of  0.5  a.u.  for  a  total 
of  9  grid  points  in  this  coordinate.  The  complete  three-dimensional  grid  contains  5049 
points.  In  Appendix  F  the  three  adiabatic  PESs  are  plotted  for  the  configurations  with  Civ 
symmetry  (9  —  n/2 )  and  Dmh  symmetry  (9  —  0). 
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This  grid  is  coarse  in  comparison  to  grids  designed  to  characterize  the  topology  of 
a  conical  intersection.  Numerical  grids  in  these  studies  are  typically  no  larger  than  1  au 
distance  from  a  conical  intersection.  '  The  coarse  grid  introduces  step-size  error  that 
enters  into  the  line  integral  used  to  calculate  the  ADT  mixing  angle.  It  also  tends  to 
under  sample  the  sharp  features  that  exist  near  conical  intersections.  An  example  of  this 
is  shown  in  Figure  1 1 . 


Figure  11.  Contour  plot  of  the  derivative  coupling  component  tr(R,  8 )  for  r  =  0.7  au  with  contours  starting 
at  0.1  au  and  spaced  in  0.195  au  increments 


In  Figure  1 1  the  contours  begin  at  0.1  au  and  are  spaced  in  increments  of  0.195  au.  The 
derivative  coupling  component  tr  rises  to  a  sharp  peak  near  the  collinear  configuration 
for  (r,  R )  =  (0.7,  8.0)  au.  However,  the  peak  appears  to  vanish  for  the  collinear 
configuration  ( 0  —  0).  A  sharp  peak  does  exist  in  the  collinear  configuration,  but  the  grid 
does  not  have  sufficient  resolution  to  capture  it. 

The  numerical  grid  for  this  work  also  has  regions  where  the  nuclei  are  in  close 
proximity.  The  electronic  structure  calculation  was  not  performed  when  two  nuclei  were 
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within  a  proximity  threshold.  The  proximity  threshold  was  set  to  0.699  au  for  H-H 
distances  and  1.17  au  for  H-B.  An  energy  ceiling  of  0.2707  au  was  established  for  these 
regions.  These  regions  have  no  impact  on  the  dynamics  of  the  system  for  small  energies. 

Aside  from  these  omissions,  the  raw  data  had  other  noticeable  blemishes  as 
shown  in  Figure  12. 


n  2. 


Figure  12.  A  slice  of  the  Vyy  adiabatic  PES  for  r  =  1.402  au  adiabatic  PES  showing  two  noticeable 
outlying  points 

Figure  12  is  a  slice  of  the  adiabatic  PES  surface  for  r  =  4.402  au.  In  this  slice  there 
are  two  obvious  outlying  points.  The  figure  also  has  a  mesa-like  structure  where  the 
proximity  threshold  of  0.2707  au  was  applied.  It  is  not  known  how  these  points 
converged  to  these  values,  nor  was  it  possible  to  recalculate  the  data  at  these  points.  The 
data  for  outlying  points  were  interpolated  using  a  low-order  polynomial  fitting  function 
and  compared  with  neighboring  values  to  ensure  continuity. 

The  raw  DCT  data  also  had  numerical  artifacts  that  were  corrected.  The  DCTs 
were  numerically  determined  to  an  arbitrary  phase  factor.  Since  the  PESs  and  DCTs  are 
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real  valued,  this  phase  factor  resulted  in  an  arbitrary  sign  being  assigned  to  a  value.  The 
correction  of  this  phase  error  is  described  by  Belcher.  An  attempt  was  made  to  identify 
outlying  points  in  the  phase-corrected  NACT  data.  This  was  abandoned  since  it  is 
difficult  to  distinguish  an  outlying  point  from  the  naturally  sharp  features  that  exist  in  the 
DCT  data.  In  Appendix  G  a  selection  of  DCT  surfaces  are  presented. 

The  Cartesian-to- Jacobi  Coordinate  Transform 

The  electronic  structure  code  required  Cartesian  coordinate  inputs.  In  Cartesian 
coordinates,  the  B  +  FE  system  has  9  DOF;  three  Cartesian  coordinates  for  each  nuclei. 
These  coordinates  can  be  transfonned  to  a  new  set  of  coordinates  which  consist  of  three 
coordinates  for  the  CM,  three  Euler  angles,  and  three  internal  coordinates.  The  internal 
DOF  of  are  described  by  the  Jacobi  coordinates  R ,  the  distance  between  the  boron  atom 
and  the  hydrogen  molecule  CM;  r,  the  bond  length  of  the  hydrogen  molecule;  and  9,  the 
polar  angle  between  the  axis  containing  the  boron  atom  and  the  axis  containing  the 
hydrogen  atoms.  This  is  shown  in  Figure  3. 

In  the  absence  of  external  fields  the  PESs  and  NACTs  are  not  dependent  on  SF 
coordinates,  and  without  loss  of  generality  the  B  +  FB  system  can  be  restricted  to  the  SF 
xz  plane.  The  orientation  of  the  B  +  FB  configuration  in  the  plane  is  given  by  the  BF 
internal  coordinates,  a  single  Euler  angle,  and  the  coordinates  of  the  B  +  FE  CM.  Table  1 
lists  the  Cartesian  and  SF/BF  coordinates  used. 
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Table  1.  The  Cartesian  and  SF/BF  coordinates  used  to  specify  nuclear  configurations  for  the  electronic 
structure  calculation 


Cartesian  Coordinates 

SF/BF  Coordinates 

XH1 

r 

Jhi 

R 

XH2 

9 

37/2 

xcm 

XB 

Ycm 

Vb 

CO 

The  relationship  between  these  coordinates  is  illustrated  in  Figure  13. 


Figure  13.  The  geometrical  relationship  between  Cartesian  coordinates  and  Jacobi  coordinates  for  the  B  + 
FC  system 


The  B  +  H2  system  lies  within  the  plane  defined  by  the  SF  X  and  Z  axes.  The  internal 
Jacobi  coordinates  r,  R,  and  6  are  the  same  described  in  Figure  3.  The  vector  rcm  is  given 
by  the  Cartesian  components  ( xcm ,  zcm)  and  points  to  the  CM  of  the  B  +  H2  system.  The 
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CM  has  been  displaced  for  clarity.  The  angle  co  is  the  polar  angle  with  respect  to  the  SF  Z 

axis.  The  transfonnation  between  these  coordinates  is  given  by  the  following  equations: 

r  =  {(xH2  ~  xH1  )2  +  (zH2  -  zH1)2}1/2 
R  —  [(A#2  +  2Xh1  —  Xb)  (iZff2  _  Zfi)  } 

a  -1  f(lXH 2  +  \XH1  ~  xb)(.XH2  ~  XH1 )  +  (|ZH2  +  |Ztfl  —  Zb)(ZH2  ~  zHl)l 

U  —  COS  i - f 

(  rR  j 


f  (|XH2  +  \XH1 

-xb)1 

1  R 

J 

1 

^  (™bxb  +  +  mHxH2) 

1 

zcm  "77  (.mBzB  "t"  ^-hzHl  "t  X^-HZH2>) 


(141) 


The  electronic  structure  calculation  is  invariant  with  respect  to  the  location  of  the 
CM  and  the  Euler  angle  co.  Hence,  only  three  coordinates  are  required  to  specify  the 
internal  configuration  of  the  B  +  H2  system.  The  numerical  calculation  was  simplified  by 
defining  a  new  set  of  coordinate  axes  shown  in  Figure  13  as  xy  and  xy.  The  boron  atom 
is  located  at  the  origin  of  the  xy  and  xy  axes.  One  of  the  hydrogen  atoms  is  fixed  on  the 
xy  axis,  while  allowing  the  other  hydrogen  atom  to  move  freely  in  the  plane  as  internal 
coordinates  changed. 

The  NACTs  also  had  to  be  transfonned  from  Cartesian  coordinates  to  Jacobi 
coordinates.  As  defined  in  Eq.  (  21  ),  the  DCTs  involve  partial  derivatives  with  respect  to 
nuclear  coordinates.  This  complicates  the  coordinate  transformation.  To  obtain  the 
NACTs  as  functions  of  internal  Jacobi  coordinates  the  Cartesian  partial  derivatives  must 
be  transfonned  to  partial  derivatives  with  respect  to  internal  Jacobi  coordinates.  An 
example  of  this  is  given  in  Eq.  (  142  )  for  the  internal  Jacobi  coordinate  R. 
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Eq.  (  142  )  uses  the  chain  rule  to  express  the  partial  derivative  operator  —  as  the  sum 


partial  derivate  operators  with  respect  to  Cartesian  coordinates.  The  chain  rule  can  be 
expressed  for  all  of  the  partial  derivatives  via  the  following  matrix  equation. 
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The  matrix  of  partial  derivatives  in  Eq.  (  143  )  transforms  the  partial  derivative  with 
respect  to  Cartesian  coordinates  to  partial  derivatives  with  respect  to  internal  Jacobi 
coordinates.  However,  the  coordinate  transformation  equations  (  Eq.  (  141  )  )  cannot  be 
inverted  analytically.  Consequently,  it  is  not  possible  to  obtain  analytic  expressions  for 
the  expressions  in  the  transfonnation  matrix  of  Eq.  (  143  ).  The  coordinate 
transformation  equations  can  be  used  to  express  the  following  chain  rule  relation. 
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The  matrices  in  Eq.  (  143  )  and  (  144  )  are  inverses  of  one  another  as  shown  in  the  next 
equation. 
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The  transformation  was  accomplished  by  using  Eq.  (  141  )  to  compute  the  matrix  of 
partial  derivatives  with  respect  to  Cartesian  coordinates.  This  matrix  was  inverted  and 
used  to  compute  the  proper  transformation  to  Jacobi  partial  derivatives. 
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Path  Dependence  of  the  APT  Mixing  Angle 

With  properly  transformed  DCTs,  the  line  integral,  Eq.  (  35  ),  can  be  carried  out 
to  determine  the  ADT  mixing  angle.  Eq.  (  35  )  was  derived  under  the  assumption  that  the 
diabatic  derivative  couplings  would  be  negligible.  This  was  enforced  by  setting  ff2  =  0 
in  Eq.  (  33  ).  The  condition  that  z^2  must  be  curl-free  (  Eq.  (  36  )  )  followed  from  this 
assumption.  Under  these  conditions  the  line  integral  is  expected  to  be  path  independent. 
However,  as  discussed  earlier,  truncating  the  electronic  basis  set  ignores  residual 
coupling  between  the  1  A ',  2  "A ',  and  1  ~A  "  electronic  eigenstates  that  span  the  P 
manifold  and  electronic  eigenstates  outside  the  P  manifold.  This  has  the  effect  of 
introducing  a  transverse  component  into  f  f2  (  Eq.  (  39  )  ). 53,56  This  component  is  not 
removed  by  the  ADT,  invalidating  the  key  assumption  used  to  derive  Eq.  (  35  ).  Hence 
the  error  introduced  by  this  component  is  termed  nonremovable  error. 

This  nonremovable  derivative  coupling  also  invalidates  the  curl-free  condition  as 
shown  by  Eq.  (  41  ).  Consequently,  the  line  integral  is  no  longer  expected  to  be  path 
independent.  To  examine  the  path  dependence  of  the  line  integral,  two  different  paths 
were  chosen.  Both  paths  use  the  same  starting  point  of  (r0,  90l  R0 )  =  (4.402, 7r/2,10) 
au.  Path  A  integrates  along  the  internal  Jacobi  coordinates  in  the  following  order:  r,  R, 
and  0.  Path  B  integrates  along  the  order  6,  r,  and  R.  Both  paths  end  at  the  point  labeled 
(r,  9,  R).  Figure  14  illustrates  these  two  paths. 
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Figure  14.  Depiction  of  Path  A  (solid  line)  and  B  (dashed  line)  within  the  3-dimensional  coordinate  space 
defined  by  the  Jacobi  coordinates  r,  0,  and  R 

Figure  15  and  Figure  16  compare  these  two  paths  for  nuclear  configurations  for  which 
9  —  n/20.  This  value  of  9  was  chosen  to  avoid  the  step  size  error  in  the  derivative 
couplings  illustrated  in  Figure  11.  Path  B  integrates  along  the  R  direction  of  for  the 
collinear  direction  and  is  influenced  by  the  step  size  error.  Figure  15  plots  the  ADT 
mixing  angle  y  as  a  function  of  the  Jacobi  coordinate  R  along  the  line  where  6  —  n/20 
and  r  =  0.7  au.  The  black  circles  in  Figure  15  indicate  the  line  integral  passed  through 
regions  where  the  electronic  structure  calculation  was  not  perfonned. 
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R  (au) 

Figure  15.  A  comparison  of  the  ADT  mixing  angle  y  calculated  from  the  line  integral  using  Path  A  (solid 
line)  and  Path  B  (dashed  line)  along  the  line  where  8  —  tt/20  and  r  =  0.7  au 


Path  A  (solid)  and  Path  B  (dashed)  are  nearly  identical  in  this  figure  indicating  that  the 
contribution  of  the  transverse  component  of  t^2  is  small  in  this  region. 

Figure  16  plots  the  ADT  mixing  angle  y  as  a  function  of  the  Jacobi  coordinate  R 
along  the  line  where  9  —  7r/20  and  r  =  2.402  au. 


Figure  16.  A  comparison  of  the  ADT  mixing  angle  y  calculated  from  the  line  integral  using  Path  A  (solid 
line)  and  Path  B  (dashed  line)  along  the  line  where  8  =  n/20  and  r  =  2.402  au 
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In  Figure  16  the  paths  are  nearly  identical  for  values  of  R  >  6.  The  paths  diverge  for 
smaller  values  of  R.  It  is  not  clear  from  the  figure  which  path  is  more  adversely  affected 
by  the  transverse  component,  nor  is  the  true  value  of  the  ADT  mixing  angle  known. 

In  general,  when  R  is  large  both  paths  are  nearly  identical  for  all  values  of  the  FF 
bond  length  r.  As  R  decreases  the  paths  diverge  as  r  becomes  large.  The  path  difference 
indicates  the  presence  of  the  transverse  component  of  Tf2-  This  was  further  confirmed  by 
calculating  the  curl  of  f f2  ■ 


0-5  e  (rad) 


Figure  17.  The  magnitude  of  the  curl  of  the  derivative  coupling  vector  field  for  r 


2.402  au 


Figure  17  shows  the  total  magnitude  of  the  curl  of  ff2  for  the  plane  where  r  =  2.402  au. 
The  curl  was  not  calculated  for  points  along  the  boundary  of  the  data  set.  The  curl  is 
non-zero  which  is  another  indication,  according  to  Eq.  (41  ),  of  the  presence  of  a  non¬ 
zero  transverse  component  of  f f2 .  However,  Eq.  (41)  does  not  indicate  the  magnitude 
of  the  transverse  component. 
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The  transverse  and  longitudinal  components  of  rf2  can  be  computed  using  the 
following  definitions  found  in  Chapter  6  of  Jackson’s  Classical  Electrodynamics  (3rd 
Ed.)68: 
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This  decomposition  was  not  performed  given  the  under-sampled  features  in  the  NACTs. 
An  accurate  determination  of  z^-f  would  require  the  electronic  structure  calculation  to  be 
perfonned  on  a  denser  numerical  grid. 
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Expected  Behavior  of  the  APT  Mixing  Angle  and  the  Diabatic  Surfaces 

The  effect  of  rf/  can  be  further  characterized  by  starting  and  ending  the  line 
integral  at  points  where  the  true  value  of  the  ADT  mixing  angle  is  known.  Using 
arguments  given  by  Alexander-  ’  ,  the  ADT  mixing  angle  can  be  predicted  for  certain 
symmetry  configurations  of  B  +  TE.  While  Alexander’s  data  was  restricted  to  the 
equilibrium  TE  bond  length  r,  these  arguments  are  valid  for  all  values  of  r. 

The  diagonal  diabatic  PESs  Vxx,  Vyy,  and  Vzz  obey  the  same  symmetry  relationship 
that  the  adiabatic  PESs  obey:  V(r,n  —  9,R)  =  V (r,  6,  R ).  The  diabatic  coupling  surface 
Vxz  is  different.  In  the  range  0  <  0  <  ^  the  center-of-charge  of  the  TE  molecule  lies 
closest  to  the  lobe  of  positive  phase  of  the  2 p  orbital  perpendicular  to  BF  Z  axis  while  for 
-  <  6  <  n  the  center-of-charge  lies  closest  to  the  lobe  of  negative  phase.  Thus,  the 
diabatic  coupling  surface  Vxz  obeys  the  symmetry  relation  V(r,n  —  9,  R )  =  —  V (r,  9,  R). 
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This  requires  that  Vxz(9  =  0)  =  Vxz(d  =  2)  =  0  to  maintain  the  continuity  o  1'  Vxz . 
According  to  Eq.  (  38  ),  this  requires  the  ADT  mixing  to  be  zero  or  n/2  for  0  =  0  and 
7t/2.  The  value  of  the  ADT  mixing  angle  is  set  by  examining  the  behavior  of  the  l  A' 
and  22 A '  PESs. 

As  stated  earlier  the  1  A '  and  2 "A '  orbitals  can  be  chosen  to  correspond  to  the  pz 
and  px  orbital  of  the  boron  atom  respectively.  The  ADT  given  by  Eq.  (31)  can  be  seen 
as  an  orthogonal  rotation  of  these  orbitals  in  the  SF  XZ  plane,  where  the  orientation  of  the 
FE  molecule  influences  the  magnitude  of  the  rotation.  When  the  B  +  FE  system  has  C2V 
symmetry  ( 9  —  n/Z)  the  1  A '  and  2  A '  adiabatic  PESs  do  not  cross.  This  is  illustrated  in 
Figure  18. 


Figure  18.  A  plot  of  the  1  2 A'  (solid  line)  and  2  2 A'  (dashed  line)  adiabatic  PESs  for  C?v  geometry  and 
r=  1.402  au 

2 

In  Figure  18  the  solid  line  corresponds  to  the  1  "A'  PES  and  the  dashed  line  corresponds 

to  the  2  A '  PES.  This  figure  can  be  compared  with  the  upper  plot  in  Figure  3  in 

26 

Alexander.  As  the  boron  atom  approaches  the  FE  molecule  in  the  perpendicular 
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2  9 

geometry  the  px  and  pz  do  not  change  orientation.  Since  the  1  A '  and  2  ~A '  adiabatic 
PESs  do  not  cross  the  1  A '  PES  maintains  the  pz  orbital  orientation  and  the  2  A ' 
maintains  the  px  orbital  orientation.  The  ADT  mixing  angle  is  constant,  here  taken  to  be 
zero  for  all  values  R  along  the  line  where  9  =  n/2  and  r  =  1.402  au. 

For  the  collinear  configuration  (9  =  0  and  r  =  1.402  au),  Cmfl  symmetry,  the  1  A' 

2 

and  2  "A '  adiabatic  PESs  appear  to  meet  near  R  =  6.5  as  shown  in  Figure  19.  Again  the 

2  2 

solid  and  dashed  lines  correspond  to  the  1  . A '  and  2  ~A '  PESs  respectively.  This  figure 
can  be  compared  with  the  lower  plot  in  Figure  3  in  Alexander.  If  the  numerical  grid  in 
Figure  19  were  more  densely  sampled  the  PESs  would  touch  at  a  single  point. 


Figure  19.  A  plot  of  the  1  2 A'  (solid  line)  and  2  2 A'  (dashed  line)  adiabatic  PESs  for  Cxh  geometry  and 
r  =  1 .402  au 

This  is  an  example  of  a  conical  intersection.  As  the  PESs  approach  R  =  6.5  au  each 
undergoes  a  sharp  change  in  slope.  The  slope  of  the  1  . A '  PES  surface  as  it  approaches 
R  =  6.5  au  from  the  left  most  closely  matches  the  slope  of  the  2  A'  PES  surface  as  it 
continues  to  greater  values  of  R.  Likewise,  the  slope  of  the  2  ~A '  PES  surface  as  it 
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2 

approaches  R  =  6.5  au  from  the  left  most  closely  matches  the  slope  of  the  1  A '  PES 

2  2 

surface  as  it  continues  to  greater  values  of  R.  The  orientations  of  the  1  A '  and  2  A '  wave 
functions  undergo  a  sharp  change  near  R  =  6.5  au.  For  R  <  —6.5  au  the  1  "A' wave 
function  corresponds  to  the p-  orbital.  For  R  <  —6.5  au  the  1  A'  wave  function 
corresponds  to  the  px  orbital.  The  2  ~A '  wave  function  undergoes  a  similar  transition. 

This  sharp  transition  is  the  source  of  the  singularity  in  the  DCTs  in  this  region. 

The  goal  of  the  ADT  is  to  remove  this  singularity.  Figure  19  suggests  that  this  can 
be  done  by  choosing  the  ADT  mixing  angle  so  that  the  resulting  diabatic  PESs  for  the 
collinear  configuration  are  associated  with  the  same  orbital.  This  eliminates  the  sharp 
change  in  orientation  thereby  eliminating  the  derivative  coupling  in  the  diabatic  basis. 
This  is  accomplished  by  setting  the  ADT  mixing  angle  equal  to  n/2  radians  where  the  px 
state  lies  energetically  below  the pz  state  (R  <  —6.5  au),  and  setting  the  angle  to  zero 
where  the  pz  state  lies  energetically  below  the  px  state  (R  >  —6.5  au).'  The  result  of  this 
transformation  is  shown  in  Figure  20. 


Figure  20.  A  plot  of  the  Vxx  (solid  line  )  and  Vzz  (dashed  line)  diabatic  PESs  for  Cmh  geometry  and 
r  =  1 .402  au 
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In  Figure  20  the  solid  and  dashed  lines  correspond  to  the  Vxx  and  Vzz  PESs  respectively. 
After  the  transfonnation  wave  functions  associated  with  Vxx  and  Vzz  correspond  entirely 
to  the  px  and  pz  orbitals  respectively.  There  is  no  sharp  orientation  change  in  the  diabatic 
wave  functions  thus  the  derivative  coupling  term  is  eliminated. 

With  the  behavior  of  the  ADT  mixing  angle  established  for  C2V  and  Cmh 
geometries  the  effect  of  x^  on  the  ADT  mixing  can  be  further  characterized  by 
perfonning  the  line  integral  along  a  path  that  begins  in  a  C2V  configuration  and  ends  in  a 
Cmh  configuration.  A  CA  configuration  is  chosen  as  a  starting  point  since  value  of  the 
ADT  mixing  angle  zero  for  all  C2V  configurations.  This  eliminates  the  requirement  of 
choosing  a  single  reference  point  for  the  line  integral.  This  path  is  illustrated  in  Figure 
21. 


Figure  21.  An  illustration  of  the  line  integral  path  chosen  to  take  advantage  of  symmetry  derived  boundary 
conditions 
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The  line  integral  is  performed  starting  at  0  —  n/2  and  ending  at  0  —  0  for  every  value  of 
R  for  a  given  value  of  the  H2  bond  length.  Figure  22  presents  the  results  of  the  ADT 
mixing  angle  calculated  using  the  path  described  in  Figure  21  for  r  =  1.402  au  and 
0  =  0.  The  dot-dash  lines  indicate  zero  and  n/2  levels  of  the  ADT  mixing  angle.  The 
black  dots  indicate  the  line  integral  passed  through  regions  where  the  electronic  structure 
calculation  was  not  performed.  Figure  22  has  a  step -function-like  transition  between  zero 
and  n/2  near  R  ~  6.5  a.u.  as  expected  from  Figure  19;  however,  the  sharp  features 
immediately  before  and  after  the  transition  are  due  to  step-size  error. 


Figure  22.  The  ADT  mixing  angle  calculated  using  the  line  integral  starting  at  8  =  7t/2  and  ending  at 
8  —  0  for  r  =  1.402  au 


Figure  22  also  shows  a  deviation  from  7r/2  for  R  >  8  au  and  a  deviation  from  zero  R  < 
3.  These  deviations  are  due  to  contributions  of  to  the  line  integral. 

Figure  23  presents  another  example  of  how  the  transverse  component  can  affect 
the  line  integral.  Here  the  line  integral  is  carried  out  for  r  =  2.402  au  and  0  =  0. 
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Figure  23.  The  ATD  mixing  angle  calculated  using  the  line  integral  starting  at  9  —  7t/2  and  ending  at 
9  =  0  for  r  =  2.402  au 

2  2 

For  this  value  of  r  and  0  the  1  A '  and  2  A '  wave  functions  do  not  experience  a  rapid 
change  in  orientation.  The  ATD  mixing  is  expected  to  remain  n/2.  As  shown  in  Figure 
23  angle  has  a  value  near  n/2;  however,  as  R  gets  smaller  it  smoothly  deviates  from 
n/2.  This  is  a  dramatic  manifestation  of  the  nonremovable  component  of  the  derivative 
coupling  vector  field. 

In  general  for  the  0  =  0  plane,  for  r  <  1.902  au  the  ADT  mixing  angle 
transitions  from  n/2  to  zero  at  some  point  as  R  becomes  smaller  for  the  collinear 
configuration.  The  jagged  features  near  these  transitions  are  due  to  step  size  error  in  the 
DCTs  as  discussed  earlier.  However,  the  smooth  deviations  from  zero  and  n/2  are  due 
to  the  nonremovable  component  of  the  derivative  coupling  error.  This  error  dominates 
the  ADT  mixing  angle  r  >  2.402  au,  and  is  illustrated  in  Figure  23. 

The  ADT  mixing  angle  passes  on  the  errors  introduced  by  f-f/to  the  diabatic 
surfaces.  This  is  clearly  seen  in  the  Vxz  diabatic  PES.  Figure  24  shows  the  Vxz  PES 
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calculated  for  r  =  2.902  au  bond  length.  The  diabatic  coupling  surface  is  not  zero  along 


the  line  9  =  0  as  predicted  by  symmetry  arguments.  Figure  24  demonstrates  that  the 
error  introduced  into  the  diabatic  PESs  can  be  energetically  significant. 


A 


0.04- 
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Figure  24.  The  diabatic  coupling  surface  Vxz  for  r  =  2.902  au  calculated  using  the  ADT  mixing  angle 
determined  by  the  line  integral  starting  at  6  =  n/2  and  ending  at  6  =  0 

The  deviation  from  zero  is  on  the  order  of  0.01  au.  The  scattering  calculation  discussed 
later  probes  energies  in  the  range  of  zero  to  0.01  au.  This  deviation  will  adversely  impact 
the  molecular  dynamics  calculated  using  these  surfaces. 

Weighted-Path  Line  Integral 

Unless  the  derivative  coupling  vector  field  z^2  is  decomposed  into  its  longitudinal 
and  transverse  components  the  contribution  of  the  transverse  component  to  the  ATD 
mixing  angle  cannot  be  removed.  As  mentioned  earlier,  this  decomposition  was  not 
attempted  in  this  work  due  to  the  under-sampled  features  in  the  NACTs.  As  the  line 
integral  is  performed  error  from  z^2  is  accumulated  along  the  path  with  the  points 
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furthest  from  the  reference  point  suffering  the  most  contamination.  The  symmetry 
derived  boundary  conditions  helped  mitigate  this  error  by  reducing  the  length  of  the  path 
required  to  make  a  determination  of  the  ADT  mixing  angle.  The  path  described  in  Figure 
2 1  involved  a  single  integration  along  a  path  in  the  r g  component  of  the  derivative 
coupling  vector  field.  This  integration  scheme,  however,  did  not  prevent  the  error 
displayed  in  Figure  24. 

The  line  integral  detennines  the  ADT  mixing  angle  up  to  an  overall  constant  y0. 

A  single  starting  point  is  selected  to  ensure  this  constant  is  the  same  for  all 
detenninations  of  the  ADT  mixing  angle  for  a  given  data  set.  The  symmetry  derived 
boundary  conditions  provide  a  way  to  begin  the  line  integral  from  multiple  points  while 
setting  the  overall  constant  consistently  across  the  data.  The  direction  of  the  integration 
in  Figure  21  can  be  reversed  provided  the  symmetry  derived  boundary  conditions  are 
used  to  set  the  constant  appropriately.  When  this  is  done  the  value  of  the  ATD  mixing 
angle  is  set  properly  in  the  9  —  0  plane  and  the  integral  accumulates  error  along  the  path 
to  the  9  =  n/2  plane.  Each  path  accumulates  the  same  magnitude  of  error  at  the  end  of 
the  path. 

For  clarity,  Path  1  is  defined  as  the  path  that  begins  in  the  9  =  n/2  plane  and 
ends  in  the  9  —  0  plane.  Path  2  begins  in  the  9  —  0  plane  and  ends  in  the  9  —  n/2  plane. 
These  two  paths  can  be  used  to  define  a  weighted-path  line  integral  using  the  following 
equation: 

<147> 
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In  Eq.  (  147  )  7i(0)  and  y2(& )  are  the  ADT  mixing  angle  calculated  along  Path  1  and  2 
respectively.  A  linear  weighting  is  used  to  emphasize  the  contribution  of  each  path  in  the 
region  where  it  has  accumulated  the  least  error.  Figure  25  illustrates  these  three  paths. 


Figure  25.  The  ADT  mixing  angle  plotted  as  function  of  6  for  Path  1  (solid  line),  Path  2  (dashed),  and  the 
weighted  path  (dash-dot)  for  r  =  2.402  au  and  R  =  3.0  au 

Figure  25  plots  the  ADT  mixing  angle  for  the  path  corresponding  to  r  =  2.402  au  and 
R  —  3.0  au.  Path  1  (solid  line)  is  set  to  the  symmetry  derived  value  of  zero  at  9  —  n/2, 
but  fails  to  reach  the  expected  value  of  y  —  n\2  at  0  —  0.  Path  2  (dashed  line)  is  set  to 
the  symmetry  derived  value  of  n/2  at  9  =  0,  but  fails  to  reach  the  expected  value  of 
y  =  0  at  6  =  n\2.  Paths  1  and  2  have  accumulated  the  same  magnitude  of  error  at  the 
end  of  each  path.  The  weighted  path  (dash-dot)  is  obtained  from  Eq.  (  147  ).  While  the 
weighted  path  has  the  correct  symmetry  derived  boundary  conditions,  it  is  not  known  to 
what  extent  the  accumulated  error  is  mitigated.  Figure  15  and  Figure  16  indicate  that  the 
transverse  component  is  not  evenly  distributed  throughout  the  data  set. 
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The  weighted-path  line  integral  approach  was  used  to  calculate  the  ADT  mixing 
angle  for  the  entire  data  set.  Only  Path  1  was  used  when  regions  where  the  electronic 
structure  calculation  was  not  perfonned  were  encountered.  Figure  26  shows  the  Vxz  PES 
calculated  for  r  =  2.902  au  bond  length  this  time  using  the  ADT  determined  by  the 
weighted-path  line  integral. 
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Figure  26.  The  diabatic  coupling  surface  Vxz  for  r  =  2.902  au  constructed  using  ADT  mixing  angle 
calculated  from  the  weighted-path  line  integral 

Unlike  Figure  24,  the  diabatic  coupling  surface  Vxz  is  now  zero  along  the  line  0  =  0.  The 
weighted-path  line  integral  successfully  corrects  the  deviations  from  the  symmetry 
derived  boundary  condition;  however,  it  is  sensitive  to  step-size  error  in  the  data.  This  is 
illustrated  in  the  next  series  of  figures.  The  diabatic  coupling  surface  Vxz  displayed  in 
Figure  27  has  a  step-size  induced  error  along  the  line  0  =  0  near  R  =  8. 
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Figure  27.  The  diabatic  coupling  surface  Vxz  for  r  =  0.7  au  constructed  using  ADT  mixing  angle  calculated 
from  the  line  integral  starting  at  9  —  n/2  and  ending  at  6  —  0 


This  diabatic  surface  was  the  result  of  an  ADT  detennined  by  integration  using  only 
Path  1.  The  derivative  coupling  surface  for  the  component  r g12  f°r  r  =  0.7  is  shown  in 
Figure  28  to  illustrate  the  source  of  this  step-size  error. 


Figure  28.  The  derivative  coupling  surface  for  the  component  xf12  for  r  =  0.7  au 


The  sharp  discontinuity  near  in  R  =  8  au  is  indicative  of  a  conical  intersection.  The  Tq  12 
component  of  the  derivative  coupling  surfaces  was  interpolated  using  a  cubic  spline  fit  in 
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an  unsuccessful  effort  to  alleviate  step-size  error.  The  grid  of  nuclear  coordinates  does 
not  sample  these  sharp  features  well  enough  for  interpolation  to  be  of  any  benefit.  The 
cubic  spline  fit  did  not  model  the  shape  of  the  discontinuity  along  6  well.  The  next  figure 
illustrates  the  effect  the  discontinuity  near  R  =  8  au  in  Figure  28  has  on  the  diabatic 
surface  constructed  from  the  ADT  determined  by  the  weighted-path  line  integral  over  the 
uninterpolated  derivative  couplings. 
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Figure  29.  The  diabatic  coupling  surface  Vxz  for  r  =  0.7  au  constructed  using  ADT  mixing  angle  calculated 
from  the  weighted-path  line  integral 

Figure  29  shows  that  the  weighted-path  line  integral  performed  well  in  the  regions  where 
the  derivative  coupling  term  is  smooth.  However,  according  to  Eq.  (  147  )  Path  2  is 
weighted  more  heavily  than  Path  1  near  0  =  0.  Path  1  is  affected  by  the  discontinuity  in 
the  DCT  and  carries  this  error  through  the  path.  All  regions  affected  by  this  step-size 
error  were  smoothed  using  cubic  spline  interpolation  as  shown  in  Figure  30. 
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Figure  30.  The  smoothed  diabatic  coupling  surface  Vxz  for  r  =  0.7  au  constructed  using  ADT  mixing  angle 
calculated  from  the  weighted-path  line  integral 

A  selection  of  diabatic  surfaces  constructed  using  the  ADT  determined  from  the 
weighted-path  line  integral  are  presented  in  Appendix  H. 

Diabatic  Potential  Energy  Surface  Fitting 

The  diabatic  PESs  are  required  to  represent  the  electrostatic  interaction  potential 
Vei  in  the  asymptotic  basis.  Eqs.  (  74  ),  (  75  ),  and  (  76  )  express  the  analytical 
representation  of  Vei  in  the  asymptotic  basis.  This  representation  includes  the  expansion 
coefficients  xan  (r>  R)  which  are  determined  from  a  numerical  fit  to  the  diabatic  PESs. 
In  Eq.  (  88  )  the  diabatic  PES  are  expressed  as  an  expansion  in  terms  of  reduced  Wigner 
rotation  matrix  elements,  Eq.  (  87  ).  The  expansion  coefficients  are  calculated 
numerically  by  fitting  each  diabatic  PES  to  the  appropriate  reduced  Wigner  rotation 
matrix.  The  relationship  between  the  expansion  coefficients  1 4rAa^(r<  R)  and  the  reduced 
Wigner  rotation  matrices  is  given  by  Eqs.  (  84  )  and  (  87  ).  The  reduced  Wigner  rotation 
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matrix  elements  are  generated  by  using  the  following  relationship  with  the  associated 


Legendre  polynomials8 


MUUKr'(f^=0) 


=  c-iy 


Ur  ~  h) 
Ur  +  /r) 


(148) 


-  PAr(cos9) 


An  analytic  expression  is  derived  for  each  expansion  coefficient.  This  process  is  shown 

A 

for  Vxz  as  an  example.  The  Vxz  PES  is  fit  using  the  (1^(6)  reduced  matrix  element  given 
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Xr  (Ar  +  1) 


Inserting  Eq.  (  149  )  into  Eq.  (  88  )  yields 


(cos9) 


(r,S,R)  =  2  ^fr.wl-Uj^IXfco*#) 


(150) 


Both  sides  of  Eq.  (  150  )  are  then  multiplied  by  Ph  (cos  0)  sin  9  and  integrated  over  9 


from  0  to  tc  to  yield 


rn 

I  Vxz(.r>  d,  R)Pl'  (cos  6)  sin  9  dO 
Jo  r 


l  v 


Pl  (cos  9)  |  >  Ph  (cos  9)  sin#  d.9 


(151) 


The  associated  Legendre  polynomials  have  the  following  orthogonality  relationship: 


2  (l  TTi)  ^ 

I  PI”'(cos0)P,7'(cos  0)  sin  0  dS  =  (2l  +  w_mySk., 


(152) 


The  orthogonality  relationship  simplifies  Eq.  (  151  )  to 
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vx7(r,R) 


Jjr+j)  fn 

[Ar(Ar  +  l)]1/2J0 


Vxz(r,  R)P/  (cos  0 )  sin  9  d6 


(153) 


Since  the  data  was  only  calculated  for  9  between  0  and  n/2  the  symmetry  properties  of 
each  surface  with  respect  to  6  —  n/2  was  used  to  extend  each  surface  to  9  —  n.  The  Vxz 
has  odd  symmetry  with  respect  to  9  —  n/2  the  values  of  Vxz  change  sign  when  reflected 
about  this  plane.  The  other  potential  energy  surfaces  have  even  symmetry  with  respect  to 
this  plane.  The  analytic  expressions  for  the  expansion  coefficients  of  these  surfaces  are 


given  by 
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2^1/2]  Vd(r,  R)Pxr(cos  9)  sin  9  d9 


These  integrals  are  performed  along  the  9  direction  for  each  (r,  R )  pair.  In  Eqs.  (  153  ) 
and  (  154  )  the  A'r  symbols  have  been  replaced  with  Ar  symbols.  This  work  calculated  fit 
coefficients  for  Ar  =  0,2,4, 6,8  to  be  directly  comparable  to  the  surfaces  calculated  by 
HIBRIDON™,  a  suite  FORTRAN  programs  used  to  calculate  Alexander’s  B  +  FF 


surfaces." 


The  number  oscillations  a  given  associated  Legendre  polynomial  experiences 
increases  as  Ar  increases.  Beyond  a  limiting  value  of  Ar  the  finite  grid  is  no  longer  able 
to  sample  the  oscillations.  This  aliasing  limits  the  order  of  the  associated  Legendre 
polynomials  that  can  be  accurately  projected  on  the  numeric  grid.  An  example  of  this  is 
shown  in  Figure  3 1 .  Although  the  solid  line  clearly  represents  P|6,  when  it  is  sampled 
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using  the  9  grid  used  in  this  work  (dashed  line)  it  appears  to  be  P4° .  The  0  grid  is  marked 
with  an  asterisk. 


Figure  31.  The  associated  Legendre  polynomial  P°6  plotted  on  a  finely  sampled  grid  (solid)  and  the  grid 
used  in  this  work  (dashed)  with  grid  points  indicated  by  an  asterisk 

Associated  Legendre  polynomials  with  orders  higher  than  Ar  —  20  are  aliased  on  the  9 
grid  used  in  this  work.  Furthennore,  the  accuracy  of  the  fit  coefficients  steadily 
decreases  as  Ar  increases  due  to  step-size  error  in  the  integration  of  Eqs.  (  153  )  and 
(  154  ).  This  is  illustrated  in  Figure  32.  In  Figure  32  the  RMS  fitting  error  for  Vzz  is 
plotted  as  a  function  of  the  fit  order  Ar  for  (r,  R )  =  (1.402, 8.0).  The  RMS  error 
decreases  for  the  first  three  values  of  Ar,  but  steadily  increases  to  the  first  alias  peak. 
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Figure  32.  RMS  fitting  error  of  Vzz  calculated  as  a  function  of  fit  order  for  (r,  R )  =  (1.402,  8.0)  au 


This  effect  can  be  reduced  by  interpolating  the  diabatic  PESs  along  9  on  a  more  densely 
sampled  grid  and  then  performing  the  fit.  Unlike  the  DCTs,  the  diabatic  PESs  vary 
slowly  along  9  as  illustrated  in  Figure  33. 


Figure  33.  The  diabatic  PES  V7Z  (asterisks)  interpolated  using  a  piecewise  cubic  Hermite  polynomial  (solid 
line)  plotted  as  a  function  of  9  for  (r,  R )  =  (1.402,  8.0)  au 
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The  asterisks  mark  the  points  where  the  electronic  structure  calculation  was  carried  out. 
The  solid  line  is  a  piecewise  cubic  Hermite  polynomial  interpolation  of  the  data.  Figure 
34  is  a  comparison  of  the  RMS  error  for  the  interpolated  and  non-interpolated  data. 
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Figure  34.  A  comparison  of  the  RMS  fitting  error  of  the  non-interpolated  data  (+)  and  the  interpolated 
data  (x) 


As  shown  in  Figure  34,  the  RMS  error  of  the  interpolated  data  (+)  is  less  affected  by  step- 
size  error  in  Eqs.  (  153  )  and  (  154  ).  The  diabatic  PESs  were  interpolated  along  the  6 
direction  before  calculating  the  fit  coefficients.  Figure  35  plots  a  selection  of  these 
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coefficients  which  is  directly  comparable  with  Figure  6  of  Alexander.” 
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Figure  35.  A  selection  of  V\rxall  expansion  coefficients  plotted  as  a  function  of  R 


For  Figure  35  the  sign  of  the  F222  coefficient  was  chosen  to  be  consistent  with 
Alexander’s  of  Vd  (  Eqs.  (  89  )  and  (  90  )  ).  The  values  of  the  expansion  coefficients  are 
comparable  to  Alexander.  The  expansion  coefficients  were  then  inserted  into  Eq.  (  92  ) 
to  yield  the  Vxr\(lll  expansion  coefficients  required  in  Eq.  (  74  )  to  represent  the 
electrostatic  interaction  potential  in  the  asymptotic  basis. 


Detennining  the  FE  Potential  Energy  Surface 

As  implied  by  Eq.  (  79  ),  the  diabatic  PESs  capture  the  effects  of  the  electrostatic 
interaction  potential  Vei  and  the  FE  PES  VH2  (r).  For  this  work  the  contribution  of  VH  (r) 
was  not  removed  from  the  adiabatic  or  diabatic  PESs.  As  a  result  the  contribution  of  the 
FE  PES  must  be  removed  from  the  V000(r,  R ),  Eq.  (  92  ),  expansion  coefficients  to 
calculate  VeL  properly.  The  FE  PES  also  appears  as  a  tenn  in  the  effective  PES  matrix 
elements  given  by  Eq.  (  97  ).  The  minimum  value  of  VHz  (r)  is  also  used  to  set  the  overall 
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offset  of  the  data  so  that  VHz  (r  =  1.402)  =  0  au.  This  allows  the  results  of  this  work  to 

26  37 

compared  with  Alexander  and  Weeks. 

The  contribution  of  VH2  (r)  can  be  isolated  by  examining  the  asymptotic  behavior 
of  the  diabtic  PESs.  The  diabatic  PESs  have  the  following  limits: 


lim  Vxx(r,R ) 

R^>  oo 


Jim  Vyy(r,  R) 


lim  Vzz(r,R) 

R—>co 


VHj 0”) 


(155) 


lim  Vxz(r,R )  =  0 

R—>co 

In  the  asymptotic  limit  the  contribution  of  Vei  goes  to  zero  leaving  only  the  contribution 
of  VH2  (r).  Consequently  the  coupling  surface  Vxz  goes  to  zero  while  the  other  diabatic 
PESs  converge  to  VHi  (r)  as  R  gets  large.  Likewise  the  fitting  function  Vs  also  coverges 
to  VH2  (r)  while  Vd  converges  to  zero.  Thus  the  expansion  coefficients  Vxrxaii  in  Eq. 

(  92  )  go  to  zero  in  the  asymptotic  limit  as  expected. 

An  accurate  determination  of  VHz  (r)  requires  the  adiabatic  PESs  to  be  calculated 
for  a  variety  of  values  of  the  FL  bond  length  in  the  asymptotic  limit  (large  R).  Figure  36 
plots  the  expansion  coefficient  V000(r,R )  nearf?  =  10  as  function  ofi?  for  r  =  1.402  au 
along  with  an  exponential  fit  to  this  data.  Figure  36  demonstrates  that  the  asymptotic 
limit  has  not  been  reached  by  R  =  10  and  that  F000  (r  =  1.402,  R)  =£  0  for  large  R.  The 
values  of  the  fit  coefficients  for  F0oo(r  —  1-402,  R )  exhibit  a  near  exponential  decay  near 


R  =  10. 
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Figure  36.  The  expansion  coefficient  Fooo (r,  R)  near  R  =  10  plotted  as  a  function  of  R  for  r  =  1.402  au  (dots) 
with  an  exponential  fit  to  these  points  (solid  line) 


The  last  eight  points  ( R  =  8.6  —  10  au)  were  fit  using  a  least  squares  fit  with  following 
exponential  fitting  function: 

/(/?)  =  ate~a2(-R~a 3-)  +  a4  (  156  ) 

The  coefficient  a4  captures  the  value  of  the  fit  coefficient  in  the  asymptotic  limit.  This  fit 

is  shown  in  Figure  36  as  a  solid  line.  The  asymptotic  values  of  F000  obtained  this  way 
are  used  to  the  global  energy  offset  of  the  PESs  so  that  VHz  (r  =  1.402)  =  0  au.  Figure 
37  compares  these  values  with  the  theoretical  Morse  potential  values  and  ab  initio 
detennined  LSTH  (Liu-Siegbahn-Truhlar-Horowitz)101, 102  values.  The  values  of  VHz  (r) 
obtained  using  the  exponential  fit  of  the  F00o  expansion  coefficient  (asterisk)  are 
comparable  to  the  Morse  potential  (dashed),  but  more  closely  follow  the  LSTH  values 
(solid). 
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Figure  37.  A  comparison  of  LSTH  (solid),  Morse  potential  (dashed),  and  the  asymptotic  expansion 
coefficient  Vo00  (asterisk)  values  for  VHz  (r) 


The  values  of  VH2  (r)  must  be  interpolated  to  match  the  grid  used  for  the  effective 
PESs  (  Eq.  (  97  )  ).  Since  VH2  (r)  for  this  work  is  only  sampled  at  nine  points,  the  method 
of  interpolation  must  be  chosen  carefully.  Figure  38  compares  three  methods  for  fitting 
the  H2  PES  in  the  asymptotic  limit  for  the  range  R  —  0.7  —  2.5  au. 


0.3 


r  (au) 

Figure  38.  A  comparison  of  three  methods  for  fitting  the  IE  PES  in  the  asymptotic  limit:  cubic  spline 
(solid),  piecewise  cubic  Flermite  (dashed),  Morse  fitting  function  (dash-dot),  F0oo  data  (asterisk) 
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Overall,  the  least  squares  fit  using  a  Morse  potential  fitting  function  (dash-dot)  had  the 
most  residual  error  since  it  was  not  constrained  to  the  data  (asterisk).  The  cubic  spline 
(solid)  and  piecewise  cubic  Hennite  polynomial  (dashed)  interpolation  curves  are  nearly 
equal  for  R  >  2.0  au.  All  three  methods  have  difficulty  fitting  the  well  region  located 
near  r  =  1.402  au.  Only  the  piecewise  cubic  Hennite  polynomial  curve  preserves  the 
value  of  zero  at  r  =  1.402  au.  The  cubic  spline  curve  shifts  the  location  of  the  minimum 
value  of  the  potential  well.  The  piecewise  cubic  Hennite  polynomial  curve  was  selected 
to  construct  VH2  (r)  in  the  asymptotic  limit.  Future  determinations  of  adiabatic  PESs  for 
other  systems  should  include  calculations  made  in  the  asymptotic  limit. 

Effective  Potential  Energy  Surfaces 

The  diabatic  effective  PESs  were  generated  using  Eq.  (  97  ).  A  FORTRAN  code 
was  developed  to  read  fit  coefficient  data;  construct  and  interpolate  the  H2  PES  VHz(r); 
interpolate  the  fit  coefficients  using  cubic  spline  interpolation  to  yield  data  for  any  user 
specified  grid  of  nuclear  coordinates;  generate  the  diabatic  effective  PESs  on  the  user 
specified  grid;  and,  diagonalize  the  diabatic  effective  PESs  matrix  at  each  value  of  (r,  R ) 
fonning  adiabatic  effective  PESs. 

A  selection  of  diabatic  effective  PESs  are  presented  in  Appendix  I.  Figure  39 
compares  the  r  =  1.402  au  values  of  the  first  14  diagonal  diabatic  effective  PESs  with 
those  produced  by  HIBRIDON™.100  Comparisons  made  with  HIBRIDON™  are  made 
using  Eq.  (  90  )  as  the  definition  of  Vd.  This  does  not  affect  the  diagonal  elements  of 
Veff  (r.  R),  nor  does  the  sign  of  the  off-diagonal  elements  affect  the  dynamics. 
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Figure  39.  The  first  14  diagonal  matrix  elements  of  Ve ?ff(r,R)  for  r  =  1.402  au:  the  results  of  this  work 
(solid  line),  the  results  produced  by  FIIBRIDON™  (dashed  line) 

The  surfaces  plotted  in  Figure  39  correspond  to  the  j  =  0,  2, 4  values  of  the  FF  rotor.  The 
potential  well  for  the  HIBRIDON™  surfaces  is  lower  in  energy;  however,  the  surfaces 
share  the  same  Ff?  rotor  energy  spectrum  as  R  becomes  large.  This  indicates  that  the 
global  offset  has  been  chosen  and  applied  properly.  The  H2  rotor  levels  are  split  by  spin- 
orbit  coupling.  Figure  40  plots  the  effective  PESs  corresponding  toy  =  0  on  a  smaller 
scale  so  that  the  differences  between  the  HIBRIDON™  surfaces  and  this  work  might  be 
more  clearly  observed. 
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Figure  40.  The  first  two  diagonal  matrix  elements  of  VRAr,  R)  for  r  =  1.402  au  corresponding  to  j  =  0: 
the  results  of  this  work  (solid  line),  the  results  produced  by  HIBRIDON™  (dashed  line) 

As  shown  in  Figure  40,  the  effective  PESs  calculated  by  HIBRIDON™  have  deeper 
potential  wells  which  occur  at  smaller  values  of  R.  This  comparative  behavior  is  typical 
of  the  effective  PESs  for  higher  j. 

Figure  41  compares  the  off-diagonal  elements  of  the  first  row  of  R)  with 

those  produced  by  HIBRIDON™. 


Figure  41.  The  off-diagonal  matrix  elements  of  Fe°/(r,  R)  corresponding  to  j  =  0  and  2  (the  single  j  —  0 
surface  is  indicated):  the  results  of  this  work  (solid  line),  results  produced  by  HIBRIDON™  (dashed  line) 
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The  single  off-diagonal  surface  corresponding  to  j  =  0  is  marked  in  Figure  41.  The  other 
surfaces  correspond  to  j  =  2.  While  the  off-diagonal  surfaces  are  similar  in  form  to  those 
calculated  by  HIBRIDON™,  there  are  significant  differences.  The  differences  observed 
in  Figure  40  and  Figure  41  will  affect  the  resulting  scattering  matrix  elements  calculated 
from  these  surfaces.  This  will  be  explored  in  the  next  section. 

Effective  PES  matrix  element  are  represented  by  an  n-by-n  matrix  in  the 
asymptotic  basis  where  n  is  the  number  of  basis  functions  (  Eq.  (  98  )  ).  This  matrix  can 
be  diagonalized  via  a  similarity  transformation  at  each  (r,  R )  value  (  Eq.  (  102  )  )  to  yield 
a  diagonal  set  of  adiabatic  effective  PESs  as  shown  below. 

VeDff(r,R)12 

U+(r,  R)  I  Vgff(r,  R)21  VeDff(r,  R)22 

VeAff(r,R)n 

o  Vgf 

The  adiabatic  effective  PESs  are  the  eigenvalues  of  the  effective  PES  matrix.  These 
eigenvalues  and  are  ordered  in  ascending  order.  A  selection  of  these  surfaces  is 
presented  in  Appendix  J. 

The  adiabatic  effective  PESs  display  avoided  crossings  similar  to  those 
encountered  in  the  electronic  adiabatic  surfaces  (Figure  19).  An  example  of  these 
crossings  is  shown  in  Figure  42. 
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Figure  42.  A  slice  along  R  —  3.6  au  of  the  first  (solid  line)  and  second  (dashed)  adiabatic  effective  PESs 

In  Figure  42  a  slice  along  R  —  3.6  au  of  the  first  and  second  adiabatic  effective  PESs  is 
plotted.  This  slice  has  surface  crossings  near  r  —  2,3,  and  3.4  au.  At  this  time  it  is  not 
known  if  these  crossings  share  similar  properties  to  the  conical  intersections  encountered 
in  electronic  adiabatic  PESs,  nor  is  it  known  how  these  crossing  influence  the  dynamics 
of  the  system. 

Scattering  Matrix  Elements 

The  diabatic  effective  PES  data  for  r  =  1.402  au  was  extracted  and  used  in 
conjunction  with  code  developed  by  Weeks  to  calculate  scattering  matrix  elements  via 
the  CPM.  The  code  uses  the  SOA  to  propagate  wave  packets  into  the  interaction  region 
and  back.  Absorbing  boundary  conditions  (BCs)  are  used  to  prevent  the  reflected  wave 
packets  from  exiting  the  propagation  grid  and  aliasing.  The  absorbing  BCs  allow  the 
smaller  propagation  grid  to  accommodate  longer  propagation  times.  Table  2  lists  the 
numerical  parameters  used  to  calculate  the  scattering  matrix  elements. 
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Table  2.  Parameters  used  in  the  ID  propagation  of  wave  packets 


Item 

Value  (au) 

Coordinate  grid  size 

256 

Coordinate  grid  spacing 

0.3 

Number  of  time  steps 

250  000 

Time  step  size 

25 

Initial  wave  packet  position 

25 

Initial  wave  packet  momentum 

±5 

Initial  wave  packet  spread 

0.5 

Amplitude  of  absorbing  BC 

0.001 

Spread  of  absorbing  BC 

15 

'in 

These  are  the  same  propagation  parameters  used  in  Weeks.  The  next  figures  compare 
scattering  matrix  elements  for  transitions  out  of  the  state  |  j,  k,ja,  o>)  =  |0, 0, 1/2, 1/2). 


Figure  43.  Transition  probabilities  as  a  function  of  energy  for  the  transition  |0,  0, 1/2, 1/2)  -> 
|0,0,l/2,l/2>:  results  from  this  work  (solid  line),  results  based  on  PESs  produced  by  HIBR1DON™ 
(dashed  line) 
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Figure  44.  Transition  probabilities  as  a  function  of  energy  for  the  transition  |0,  0, 1/2, 1/2)  -> 
|0,0,l/2,3/2>:  results  from  this  work  (solid  line),  results  based  on  PESs  produced  by  HIBRIDON™ 
(dashed  line) 


The  scattering  matrix  elements  shown  Figure  43  and  Figure  44  are  characterized  by  broad 
smoothly  varying  Stueckelberg  oscillations  which  are  interrupted  by  rapidly  oscillating 
Feshbach  resonances104  near  energy  values  of  0.0015  and  0.0055  au.  While  the  transition 
probabilities  have  similar  structure  qualitatively,  these  figures  show  that  both  the 
Stueckelberg  oscillations  and  the  Feschbach  resonances  are  sensitive  to  changes  in  the 
input  diabatic  effective  PESs. 

The  Feshbach  resonances  occur  at  energies  corresponding  to  the  vibrational 
eigenvalues  of  the  weakly  bound  B  -"Hi  van  der  Waals  complex.  As  shown  in  Figure 
45  the  Feshbach  resonances  are  shifted  to  higher  energies  than  those  based  on 
HIBRIDON™  surfaces.  This  is  consistent  with  the  difference  in  potential  well  depths  of 
the  diagonal  effective  PESs.  On  average  the  potential  well  depth  for  the  first  20  diagonal 
diabatic  effective  PESs  in  this  work  is  2.6  X  10-5  au  higher  in  energy  than  the 
HIBRIDON™  data  set.  This  has  the  effect  of  raising  the  vibrational  eigenvalues  of  the 
weakly  bound  B  •••Ha  van  der  Waals  complex. 
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Figure  45.  A  plot  of  the  Feshbach  resonance  region  of  the  transition  the  transition  |0, 0, 1/2, 1/2)  -> 
|0,0,l/2,l/2):  results  from  this  work  (solid  line),  results  based  on  PESs  produced  by  HIBRIDON™ 
(dashed  line) 

The  shape  of  the  Feshbach  resonances  is  also  different.  This  is  largely  due  in  to  the 
differences  in  the  off-diagonal  coupling  terms  (Figure  41).  The  Feshbach  resonances  are 
also  sharper  for  the  HIBRIDON™  data  set.  This  suggests  that  the  shallower  well 
structure  observed  in  Figure  40  does  not  trap  evolving  wave  packet  amplitude  for  as  long. 
The  wave  packet  amplitude  exiting  the  interaction  region  makes  contributions  to  the 
correlation  function  (see  Eqs.  (  124  )  and  (  138  ))  for  smaller  durations  leading  to 
broader  features  in  the  scattering  matrix  elements. 

The  shapes  of  the  potential  wells  influence  the  phase  of  the  Stueckelberg 
oscillations.  As  seen  in  Figure  43  and  Figure  44,  the  maxima  of  the  Stuekelberg 
oscillations  calculated  in  this  work  occur  at  larger  energies  compared  to  those  of  the 
HIBRIDON™  data  set.  As  the  reactant  wave  packet  enters  the  interaction  region  it 
experiences  an  effective  PES  with  a  different  shape  than  the  surfaces  based  on  the 
HIBRIDON™  data  set.  During  the  interaction  the  evolving  wave  packet  explores 
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different  regions  of  the  other  effective  PESs  before  exiting  the  interaction  region.  The 
differences  in  the  shape  of  the  effective  PESs  is  primarily  responsible  for  the  phase 
differences  in  the  Stuekelberg  oscillations.  The  complicated  interaction  of  the  evolving 
wave  packet  with  multiple  effective  PESs  makes  it  difficult  to  predict  a  priori  how 
changes  in  the  effective  PESs  will  affect  the  resulting  scattering  matrix  elements.  It  is 
also  not  known  how  changes  in  the  effective  PESs  affect  the  calculated  cross  sections  and 
reaction  rates. 
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IV.  Conclusion 


The  diabatic  and  adiabatic  effective  PESs  for  the  inelastic  collision  of  a  boron 

2  2 

atom  with  a  hydrogen  molecule  (  Eq.  (  1  )  )  were  calculated  from  the  1  A',  2  . A and  A" 
adiabatic  electronic  PESs  and  the  associated  DCTs.  The  processing  steps  leading  to  the 
effective  PESs  are  outlined  in  Table  3. 


Table  3.  Processing  steps  for  obtaining  effective  PESs 

_ Processing  Steps _ 

1.  Ab  initio  calculation  of  electronic  adiabatic  PESs  and 
DCTs 

2.  Smooth  electronic  adiabatic  PESs 

3.  Apply  Cartesian-to-Jacobi  coordinate  transformation  to 
derivative  couplings 

4.  Calculate  ADT  mixing  angle  using  weighted-path  line 
integral 

5.  Construct  electronic  diabatic  PESs 

6.  Smooth  step-size  induced  error  in  electronic  diabatic  PESs 

7.  Fit  electronic  diabatic  PESs  to  associated  Legendre 
polynomials 

8.  Isolate  asymptotic  behavior  and  set  global  data  offset 

9.  Interpolate  fit  coefficients  to  user  defined  grid 

10.  Construct  diabatic  effective  PESs 

1 1 .  Diagonalize  diabatic  effective  PES  matrix  elements  to 
construct  adiabatic  effective  PESs 


The  adiabatic  electronic  surfaces  were  calculated  at  the  SA-MCSCF/CI  level  including 
the  H2  bond  length  as  a  DOF.29' 45-46  The  derivative  couplings  were  calculated  directly 
from  the  SA-MCSCF/CI  wave  functions  using  analytic  gradient  techniques.29,47  The 
electronic  structure  calculation  was  not  attempted  for  nuclear  configurations  with 
internuclear  distances  below  a  set  threshold.  In  addition  to  these  regions,  the  raw  data  set 
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contained  points  which  did  not  converge  to  values  near  their  neighboring  points.  In  the 
process  of  perfonning  the  electronic  structure  calculation  a  metric  based  on  the  continuity 
of  the  angular  momentum  of  the  system  was  used  to  gauge  convergence.  There  were  a 
small  number  of  points  that  met  this  metric  but  did  not  converge  properly.  These 
outlying  points  were  smoothed  in  the  electronic  adiabatic  surfaces;  however,  no  attempt 
was  made  to  smooth  the  derivative  coupling  surfaces  due  to  their  naturally  spiky 
topology.  Each  outlying  point  was  identified  individually  and  smoothed. 

After  the  DCTs  were  transformed  from  Cartesian  to  Jacobi  coordinates  the  ADT 
mixing  angle  was  calculated  by  perfonning  a  line  integral  through  the  derivative  coupling 
vector  field  (  Eq.  (  35  )  ).  When  the  ADT  mixing  angle  was  first  calculated  the  standard 
approach  of  choosing  a  single  reference  point  for  the  line  integral  was  used.  The 
resulting  line  integral  was  not  path  independent  as  assumed  when  deriving  Eq.  (  35  ). 
Residual  coupling  between  boron  2P  states  and  other  states  not  considered  in  this  work 
give  rise  to  a  transverse  component  of  the  derivative  coupling  field.  ’  This  component 
introduces  an  error  which  accumulates  along  the  line  integral  path — the  nonremovable 
error. 

This  error  has  been  characterized  in  regions  in  close  proximity  to  conical 
intersections  by  integrating  over  a  closed  path.  '  When  the  path  does  not  enclose  a 
conical  intersection  the  path  integral  should  be  zero.  When  a  conical  intersection  is 
enclosed  the  resulting  mixing  angle  should  be  nn  where  the  integer  n  is  the  number  of 
enclosed  conical  intersections.  ’  The  deviation  from  these  values  is  a  measure  of  the 
nonremovable  error. 
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This  work  represents  the  first  effort  to  characterize  this  error  for  a  wide  range  of 
nuclear  coordinates.  The  path  independence  of  the  line  integral  is  one  measure  of  the 
influence  of  the  nonremovable  error.  As  shown  in  Figure  15,  not  all  regions  in  the  data 
are  strongly  affected  by  this  error.  Previously  published  results  on  B  +  H2  for  r  =  1.402 
au  lie  in  a  region  weakly  affected  by  this  error  (Figure  22).'  ’ '  Path  integrals  in  regions 
where  r  >  1.902  au  have  significant  path  dependence  (Figure  16).  The  curl  of  the 
derivative  coupling  vector  field  further  verifies  the  presence  of  a  non-zero  transverse 
vector  field.  However,  the  path  independence  of  the  line  integral  is  not  an  absolute 
measure  of  the  nonremovable  error.  In  this  work,  the  symmetry  of  B  +  H2  is  used  to 
predict  the  true  value  of  the  ADT  mixing  angle,  providing  an  absolute  measure  of  the 
contribution  of  the  nonremovable  error  for  the  collinear  and  perpendicular  nuclear 
configurations. 

Symmetry  derived  boundary  conditions  allow  the  line  integral  to  begin  at  any 
point  at  which  the  ADT  mixing  angle  is  known.  This  allows  the  overall  length  of  the 
path  to  be  reduced  thereby  limiting  the  amount  of  accumulated  error.  The  ADT  mixing 
angle  was  calculated  by  integrating  over  the  6  component  of  the  derivative  coupling  field 
and  using  boundary  conditions  to  set  the  proper  offset.  This  technique  does  not  prevent 
significant  error  from  accumulating  (Figure  23).  This  error  is  passed  on  to  the  diabatic 
PESs  (Figure  24). 

The  ADT  mixing  angle  can  only  be  predicted  for  specific  nuclear  configurations 
with  special  symmetry.  Furthennore,  the  nonremovable  error  is  not  equally  distributed 
through  the  data  set  (Figure  17).  This  prevents  the  contribution  of  the  nonremovable 
error  from  being  removed  at  each  point  along  the  path.  For  the  B  +  H2  system,  the 
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symmetry  derived  boundary  conditions  allow  two  paths  to  be  taken  through  the  same 
data.  The  results  of  these  line  integrals  can  be  averaged  giving  weight  to  regions  less 
affected  by  the  nonremovable  error.  While  this  does  not  eliminate  the  contribution  of  this 
error,  the  weighted-path  line  integral  enforces  the  proper  boundary  conditions  (Figure 
26).  This  work  introduces  the  weighted-path  line  integral  as  a  method  of  producing  more 
accurate  diabatic  PESs.  The  technique  shows  promise  for  systems  with  a  high  degree  of 
symmetry  and  can  aid  the  calculation  of  diabatic  PESs  for  systems  with  a  larger  number 
of  nuclei. 

The  weighted  line  integral  is  also  sensitive  to  step-size  error  introduced  by  the 
finite  grid  of  coordinates.  The  derivative  coupling  surfaces  are  most  affected  by  this 
error  due  to  the  sharp  features  that  they  can  possess  (Figure  28).  The  nuclear  grid  used  in 
this  work  under  sampled  these  sharp  features  (Figure  1 1).  Consequently,  efforts  to 
reduce  this  error  by  interpolating  the  derivative  coupling  data  set  were  unsuccessful.  Step 
size  error  was  smoothed  by  masking  the  points  and  using  cubic  spline  interpolation. 

The  smoothed  diabatic  PESs  were  fit  using  associated  Legendre  polynomials 
which  in  turn  yielded  the  expansion  coefficients  required  to  construct  the  matrix  elements 
of  the  electrostatic  interaction  potential  Kf  The  asymptotic  behavior  of  the  system  was 
modeled  by  using  a  least  squares  fit  with  an  exponential  fitting  function  for  R  =  8.6  —  10 
au.  The  predicted  FE  PES  surface  agrees  with  the  LSTH  surface.101, 102  Piecewise  cubic 
Hennite  polynomial  was  used  to  interpolate  the  FE  PES  in  the  asymptotic  limit.  The 
resulting  surface  was  used  to  account  for  the  contribution  of  the  hydrogen  molecule  in  the 
full  Hamiltonian  (  Eq.  (  60  )  )  as  well  as  subtract  its  influence  in  the  F0oo  expansion 
coefficient  (  Eq.  (  92  )  ). 
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A  FORTRAN  code  was  developed  to  interpolate  the  expansion  coefficients  to  a 
user  defined  grid.  The  diabatic  effective  PESs  (  Eq.  (  97  )  )  were  calculated  on  this  grid. 
20  basis  functions  (all  basis  up  to  and  including  j  =  6)  were  used  to  compute  the  matrix 
elements  V‘t ?ff(r,  /7)^'<y  The  data  corresponding  to  r  =  1.402  au  was  extracted  from 
these  matrix  elements  and  compared  to  values  calculated  by  HIBRIDON™.  The  surfaces 
showed  good  agreement  in  shape  and  magnitude.  The  matrix  elements  Re^( r,  7?)^ 
were  diagonalized  at  each  (r,  R )  value  to  yield  adiabatic  effective  PESs.  These  surfaces 
exhibit  surface  crossings  similar  to  those  encountered  in  the  electronic  adiabatic  surfaces 
(compare  Figure  19  with  Figure  42);  however,  it  is  not  known  if  they  share  similar 
properties  associated  with  conical  intersections. 

Finally,  the  values  of  the  diabatic  effective  PESs  for  r  =  1.402  au  were  used  in 
the  same  one  dimensional  CPM  method  code  developed  by  Weeks  et  al.  to  compute 
scattering  matrix  elements.  The  resulting  scattering  matrix  elements  exhibit  Stuekelberg 
oscillations103  and  Feshbach  resonances104  shifted  to  higher  energies  than  those  based  on 
the  HIBRIDON™  data  set.  The  scattering  matrix  are  sensitive  to  the  form  of  the  input 
effective  PESs  (Figure  43  and  Figure  44).  The  difference  in  the  well  depths  of  the 
diabatic  effective  PESs  affect  both  the  Stuekelberg  oscillations  and  Feshbach  resonances. 

Summary  of  Key  Contributions 

This  work  presents  the  adiabatic  PESs  and  DCTs  of  the  B  +  H2  system  calculated 
over  an  extended  range  of  nuclear  coordinates  over  all  three  internal  nuclear  DOF  (boron 
distance,  H2  orientation,  and  H2  bond  length).  Previous  ab  initio  calculations  of  the 
adiabatic  PESs  for  the  B  +  H2  system  fixed  the  molecular  hydrogen  bond  length  at  the 
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equilibrium  value  of  1 .402  a.u.  26, 11 ' 36, 37  The  DCTs  were  used  to  calculate  the  ADT 

26  27 

mixing  angle.  Previous  used  other  methods  to  estimate  the  ADT.  ’  The  error 
introduced  by  employing  a  truncated  set  of  adiabatic  states  was  characterized  by 
examining  the  path  dependence  of  the  line  integral  and  using  symmetry  derived  boundary 
conditions.  The  weighted-path  line  integral  was  employed  as  a  new  method  for  reducing 
the  effect  this  error  has  the  resulting  diabatic  PESs.  The  procedure  for  determining  the 
Vxr\  n  expansion  coefficients  was  modified  to  account  for  variable  TE  bond  length.  The 
two  dimensional  diabatic  and  adaiabtic  effective  PESs  were  calculated.  The  effective 
PESs  data  was  extracted  for  r  =  1.402  au  and  used  to  calculate  scattering  matrix 
elements  using  the  CPM.  These  results  were  compared  with  previous  results  in  the  first 
attempt  to  observe  the  sensitivity  of  this  calculation  to  the  input  electronic  structure  data. 


Recommendations  for  Future  Work 

This  work  lays  the  foundation  for  many  significant  research  endeavors.  Table  4 
summarizes  the  significant  extensions  of  this  work. 


Table  4.  A  list  of  further  research  topics 

_ Future  Work _ 

1 .  Relax  centrifugal  sudden  approximation 

2.  Relax  pure  precession  approximation 

3.  Extend  asymptotic  basis  beyond /  =  1/2 

4.  Implement  2D  propagation  in  CPM 

5.  Examine  the  sensitivity  of  scattering  matrix  elements, 
cross  sections,  and  reaction  rates  to  input  PESs 

6.  Allow  vibrational  transitions  (v  =£  v') 

7.  Calculate  diabatic  B  +  H2  surface  for  r  >  1.902  au  using 
non-NACT  techniques 


114 


The  CS  approximation  was  implemented  in  this  work  to  reduce  the  computational 
overhead.  Relaxing  the  CS  approximation  will  require  significant  changes  to  the  code 
used  to  compute  the  effective  PESs  and  calculate  the  scattering  matrix  elements.  Once 
implemented,  a  comparison  can  be  made  between  CS  and  non-CS  scattering  matrix 
elements  to  determine  the  error  introduced  by  this  approximation. 

Relaxing  the  pure  precession  approximation  (constant  spin-orbit  coupling 
coefficient  <f)  requires  a  new  set  of  electronic  adiabatic  PESs  to  be  calculated.  This  effect 
must  be  included  as  a  relativistic  correction  in  the  electronic  structure  code.  Once 
accomplished,  the  error  introduced  by  the  pure  precession  approximation  can  be  assessed. 

The  equations  for  cross-sections  and  reaction  rates  (  Eqn.  (  139  )  and  (  140  )  ) 
require  the  asymptotic  basis  to  be  extended  beyond  a  total  angular  moment  of  J  —  1/2. 
Although  the  asymptotic  Hamiltonian  is  diagonal  with  respect  to  total  angular 
momentum,  extending  the  basis  set  will  greatly  increase  the  computational  resources 
required  to  compute  each  total  angular  momentum  block. 

The  computational  requirements  will  also  greatly  increase  when  implementing  2D 
propagation  in  the  CPM  code.  Once  implemented  the  dynamical  effects  of  the  surface 
crossings  observed  in  Figure  42  can  be  studied.  There  is  further  interest  in  observing 
how  the  added  DOF,  r,  affects  the  Feshbach  resonances  and  Stuekelberg  oscillations 
observed  in  the  transition  probabilities  (Figure  43  and  Figure  44).  2D  propagation  could 
also  serve  as  another  measure  of  the  sensitivity  of  the  scattering  matrix  calculation.  It 
would  also  be  possible  to  study  the  sensitivity  of  the  scattering  matrix  calculation  by 
examining  an  analytic  potential,  such  as  the  Morse  potential.  The  parameters  of  an 
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analytic  potential  could  be  varied  and  the  effect  on  the  resulting  scattering  matrix  could 
be  studied. 

Once  the  2D  propagation  code  has  been  implemented,  vibrational  transitions  can 
be  considered.  Considering  energies  in  this  range  also  opens  up  the  possibility  that  boron 
will  react  forming  BH2.  The  full  treatment  of  the  reaction  requires  PESs  of  the  BH  +  H 
complex;  however,  reaction  probabilities  can  be  estimated  by  measuring  the  probability 
amplitude  that  exits  the  propagation  grid  for  small  R . 

The  results  of  this  work  indicate  that  the  nonremovable  error  made  a  greater 
contribution  to  the  ADT  mixing  angle  for  H2  bond  lengths  r  >  1.902  au.  The  B  +  TE 
surfaces  calculated  by  Alexander  were  calculated  for  r  >  1.902  au.  ’  “  Thus  the 
calculation  of  the  B  +  TE  diabatic  PESs  for  r  >  1.902  au  using  methods  that  do  not  use 
DCTs  would  be  valuable  for  determining  how  these  methods  are  affected  by  the 
nonremovable  error.  This  would  also  provide  a  way  to  further  judge  the  utility  of  the 
weighted-path  line  integral  approach  to  calculating  the  ADT  mixing  angle. 
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Appendix  A  The  Antihermitian  Property  of  the  Derivative  Coupling  Terms 

The  DCTs  given  by  Eq.  (  21  )  can  be  shown  to  be  antihermitian.  This  important 
property  reduces  the  number  of  DCTs  that  must  be  computed.  The  derivation  of  this 
property  begins  by  taking  the  gradient  with  respect  to  nuclear  coordinates  on  both  sides 
of  Eq.  (  14  ),  the  orthononnality  condition  of  the  electronic  eigenstates. 

V<,n  {J  dqe<P*(qn,  qe)<Pj(qn,qe) }  = 

The  gradient  V(?ji  will  pass  through  the  integral  and  after  using  the  product  rule  the 
following  expression  is  obtained: 

j  dqe{Vqnd>*(qn,qe)}^j(qn,qe)  +  J  dqe&* (qn,  qe){Vqn<Pj(qn,  qe)}  =  0 

When  Eq.  (  21  ),  the  definition  of  the  derivative  coupling  tenn  f-j,  is  substituted  into 
this  expression  the  antihennitian  relationship  is  obtained: 
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Appendix  B  The  Derivation  of  the  Generalized  Hellmann-Feynmann  Theorem 


The  generalized  Hellmann-Feynmann  theorem,  Eq.  (  27  ),  gives  the  clearest 
indication  when  a  system  will  behave  nonadiabatically.  It  shows  that  when  the  difference 
between  the  adiabatic  PESs  for  two  different  electronic  eigenstates  is  small,  the 
derivative  coupling  f  fj  will  be  large,  possibly  singular  if  the  PESs  are  equal. 

The  Hellman-Feynman  theorem  establishes  the  following  relationship 


d 
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where  H(A)  is  a  Hermitian  operator,  \4>U))  is  a  normalized  eigenstate,  and  £’(A)  is  the 
corresponding  eigenvalue  all  of  which  depend  on  the  real  valued  parameter  A.106  This 


equation  involves  a  single  eigenstate  and  eigenvalue,  thus  if 


(W)\±H(A)\HA)) 


were 


expressed  as  a  matrix  it  would  be  diagonal.  To  compare  two  different  eigenstates  and 
eigenvalues  the  generalized  Hellmann-Feynman  must  be  used. 

The  electronic  TISE  defined  by  Eq.  (  12  )  serves  as  the  starting  point  for  this 
derivation.  In  the  expression  above,  A  is  a  collection  of  real  parameters.  These 
parameters  correspond  to  the  nuclear  coordinates  qn  in  Eq.  (  12  ).  The  derivative 

operator  ^  corresponds  to  VQn.  The  derivation  proceeds  as  follows: 
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Heleclj (*7n))  —  (*7n)  L/(*?n)) 

(i(9n)|^tec|7(<?n))  =  (iCqfJl^CqJlKflfn)) 

=  VjA(,qn)(i(,qn)\jXqn)) 

=  Vfiqn)^ 

The  derivatives  with  respect  to  nuclear  coordinates  are  then  taken  on  both  sides  and  the 
product  rule  is  used  to  expand  the  left  hand  side  of  the  equation. 

—  \y Qn^^ri)  l}{^eiec \j (flfn))}  "h  {^(.Qn)\^qn^elec\j '(flfn)) 

+  {<i(<?n)l^etec}{VqnL/(qfn))} 

The  definition  of  the  DCTs  fd,  Eq.  (  21  ),  is  then  used  to  further  simplify  the  expression. 

(KRn)\VqnBelec\j(.Rn))  +  +  ^(<7n)Ay  =  0 

After  using  antihermitian  property  of  the  DCTs  derived  in  Appendix  A,  the  general 
Hellman-Feynman  theorem  is  given  by 

->4  _  (t(flfn)  |^qn^eJec  [/ (flfn)) 

^  “  r/Oln)  -  ^(qn) 
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Appendix  C  The  Effects  of  Cs  Symmetry  on  Coupling  Terms 


Eq.  (  29  )  uses  the  result  that  only  orbitals  with  the  same  symmetry  will  have  non¬ 
zero  coupling  tenns.  This  can  be  understood  by  examining  the  fonn  of  the  coupling 
terms  within  the  context  of  Cs  symmetry.  As  shown  in  Figure  1  and  Figure  2,  functions 
with  A '  Cs  are  even  functions  with  respect  to  inversion  of  electronic  coordinates  while 
functions  A  "  Cs  are  odd  functions  with  respect  to  the  same  inversion.  An  integral  of  an 
odd  (A  "  Cs)  function  over  a  symmetric  interval  about  the  origin  is  zero.  In  Eqs. 

(  21  )  and  (  22  )  both  integrals  are  over  electronic  coordinates;  however,  the  operators  VQn 
and  Vgn  involve  derivatives  with  respect  to  nuclear  coordinates.  If  the  resulting  function 
under  the  integral  sign  is  odd,  then  the  integral  will  be  zero. 

It  can  be  shown  that  the  operators  V£?n  and  do  not  change  the  symmetry  of  a 
function.  An  example  of  this  is  expressed  in  the  following  equation: 

(*7n)  —  /A'(«In) 

The  expression  states  that  a  function  with  A 1  operated  on  by  VQn  yields  a  function  with  A 1 
symmetry.  If  (hA'( qn )  were  expanded  in  a  Taylor  expansion  each  tenn  of  the  Taylor 
expansion  must  have  symmetry  properties  consistent  with  the  symmetry  of  the  entire 
function.  A  simple  example  of  this  is  the  sine  and  cosine  functions.  A  Taylor  expansion 
of  the  sine  function  yields  only  tenns  with  odd  symmetry  while  an  expansion  of  the 
cosine  functions  yields  only  terms  with  even  symmetry.  In  the  case  of  electronic 
eigenstates  like  dv(qn),  a  multivariable  Taylor  expansion  will  result  in  various  products 
of  nuclear  and  electronic  coordinates.  When  the  operators  VQn  and  V2qn  are  applied  to  the 
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expansion,  the  electronic  coordinates  will  be  treated  as  constants  and  thus  will  not  be 
affected  by  the  operation.  This  means  the  symmetry  of  each  tenn  in  the  expansion  will 
remain  the  same  with  respect  to  electronic  coordinates.  Thus  the  function  resulting  in  the 
operation  maintains  the  same  symmetry  properties  with  respect  to  electronic  coordinates 
as  the  original  eigenstate. 

Given  that  the  functions  resulting  from  the  operators  Vq  and  Vqn  maintain  the 
same  symmetry  as  the  original  function,  it  is  straightforward  to  determine  symmetry  of 
the  function  being  integrated  in  Eqs.  (  21  )  and  (  22  ).  The  product  of  two  functions 
which  both  have  either  A '  Cs  or  A  "  Cs  will  result  in  a  function  with  A '  Cs  symmetry  and 
the  resulting  integral  will  be  non-zero.  The  product  of  a  function  with  A '  Cs  and  a 
function  with  A  "  Cs  will  result  in  a  function  with  A  "  Cs  and  the  resulting  integral  will  be 
zero.  Thus  coupling  terms  involving  eigenstates  with  different  symmetry  will  be  zero. 
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Appendix  D  The  Derivation  of  the  ADT  Mixing  Angle 


The  ADT  transformation  given  by  Eq.  (31)  introduced  the  ADT  mixing  angle 
y(qfn).  The  mixing  angle  must  be  chosen  so  that  the  derivative  couplings  in  the  diabatic 
basis  f fj  are  negligible.  The  definition  of  derivative  couplings  in  the  adiabatic  basis  f  d 
are  given  by  Eq.  (  2 1  )  as 

?ij  =  J  dQe^f*(Qn>Qe)%n^f(Qn,Qe) 

The  functions  <t>f*(qfn,  qe)  and  4>d(qrn,  qe )  will  be  chosen  to  be  real-valued  functions. 
By  choosing  real-valued  functions,  the  antihennitian  property  of  the  derivative  couplings 
(  Eq.  (  24  )  )  implies  that  diagonal  DCTs  are  zero 

f " T-  0 

•'ll  "'ll  u 

Furthennore,  the  states  with  A  'Cs  symmetry  will  not  mix  with  the  state  with  A  "Cs 
symmetry.  This  eliminates  all  but  the  t^2  and  t21  DCTs  which  are  related  to  each  other 
via  the  antihennitian  property. 

First,  the  transfonnation  given  in  Eq.  (  3 1  )  is  inverted  to  yield  the  adiabatic  states 
in  terms  of  the  diabatic  states. 

(flfni  Re)  =  (qn,  qe)  cos  Y (qn)  +  ^22a'  Re)  sin  Y (<?n) 

®22A'(qn’qe)  =  -^A'^Rrv  Re)  Siny(_qn)  +  <t>22A'(Rn’Re)C0SY(qn) 

To  allow  for  a  more  compact  notation  the  functional  dependence  on  nuclear  and 
electronic  coordinates  is  assumed.  These  equations  are  inserted  into  Eq.  (  2 1  )  to  yield 
the  following  integral 
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Ti2  =  j  dr{ <fr°2A'  C0SYiRn)  +  ®22a'  sin^(<?n)}V£?n{-^^  siny(qn) 

+  <S>2 2a'  cosy(qn)} 

When  multiplied  out,  the  integral  has  eight  terms  as  shown  below 
T12  =  J  dr[-<P ?2Ai^qn^2Ai  sin  y(q„)  cosy(qn)  -  (<bf2yl,)2  cos 2y(qn)  (v^yCqj) 

+  V^,)  cos2  y(qj 

-  siny(qfn)  cosy(qn)  (v„ny(qn)) 

-  sin2  y(qn) 

-  d^/d^/  siny(qn)  cos  y(qn)  (v^yCqj) 

+  d^W^cb^,  sin  y(qn)  cosy(qn) 

-  (d5^^)2  cos2  y(qn)  (v^yCqj)} 

The  tenns  that  depend  only  on  nuclear  coordinates  can  be  pulled  outside  the  integral. 

The  terms  Cb^W^d*^/  and  d*^' V^d^y  lead  to  diagonal  DCTs  when  integrated 
and  are  thus  zero.  When  the  orthonormailty  condition,  Eq.  (  14  ),  is  applied  tenns  like 

(^iw)  w'd  integrate  to  unity  while  terms  like  d’^/d’^/  will  integrate  to  zero.  When 
these  simplifications  are  made  the  expression  above  becomes 

T12  =  -  cos2  y(q„)  (v^yCqj)  +  cos2  y(qn)  ff2  -  sin2  y(qn)  ?£. 

-  sin2  y(qn)  (v^yCqj) 
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After  using  the  antihermitian  property  of  the  DCTs  and  applying  standard  trigonometric 
identities  the  expression  simplifies  to 

T12  =  -  cos2  y(qn)  =  -Vqny(qn)  +  ff2 
Thus  when  ff2  is  set  to  zero  Eq.  (  34  )  is  obtained. 

VrCqn)  =  -T12 

The  equation  gives  rise  to  the  line  integral,  Eq.  (  35  ),  yielding  the  ADT  mixing  angle 
y(qn)  for  each  nuclear  configuration. 
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Appendix  E  The  Symmetry  of  the  Electrostatic  Interaction  Potential 


The  symmetry  of  the  FE  molecule  influences  the  structure  of  the  electrostatic 
interaction  potential  when  represented  in  the  asymptotic  basis.  Because  H2  is  a 
homonuclear  diatomic  molecule,  when  considering  H2  in  the  BF  frame  with  the  Jacobi 
coordinate  9,  all  functions  containing  this  variable  must  have  the  property  f(n  —  0)  — 
/(0).  Thus  when  considering  the  diabatic  PESs  Vxx  ,  Vyy,  and  Vzz  given  by  Eqs.  (  88  ) 
only  expansion  coefficients  V^aiu(r,  R)  where  Xr  is  even  will  be  non-zero.  The 

coefficients  for  odd  values  of  Xr  vanish  due  to  the  following  propery  of  the  reduced 

26 

Wigner  rotation  matrix  elements”  : 

-  9)  = 

For  Vxx  ,  Vyy,  and  Vzz,  fi  is  either  0  or  1 .  Thus  when  Xr  is  odd  the  sum  /r  +  Xr  is  also  odd 
causing  a  change  in  polarity. 

Similarly,  the  coefficients  for  odd  values  of  Xr  for  Vxy  will  also  vanish  due  to  the 
requirement  that  Vxy  vanish  when  9  =  n/2.  This  requirement  is  a  consquence  of  the 
polarity  of  the  boron  2 p  orbitals  coupled  by  Vxy.  As  described  in  detail  by  Alexander  , 
the  sign  of  Vxy  will  flip  as  6  ->  n  —  9.  Given  that  /x  is  unity  for  Vxy,  odd  values  of  Xr  lead 
to  an  even  sum  /j.  +  Xr.  This  does  not  change  the  polarity  as  required.  Therefore,  all 
expansion  coefficients  for  odd  Xr  vanish  for  Vxy  as  well. 

Knowing  that  Xr  will  on  take  on  even  values,  the  property  of  3 -j  symbols 
described  by  Eq.  (81)  will  determine  which  angular  momentum  states  j  couple  with  one 
another.  Since  Xr  will  always  be  even,  /  and  j  must  both  be  either  odd  or  even  for  to 
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have  a  non-zero  expansion  coefficient.  As  a  result,  there  are  two  subsets  of  states:  even 
valued  j  states,  and  odd  valued  j  states.  States  from  one  of  these  subsets  will  not  couple 
with  a  state  in  the  other  subset.  Therefore,  states  with  even  and  odd  j  can  be  considered 
separately.91  States  with  even  values  of  j  correspond  to  parahydrogen  and  odd  values 
correspond  to  orthohydrogen.  For  this  work  only  parahydrogen  states  are  considered. 
This  reduces  the  number  of  states  that  can  be  included  in  the  basis  due  to  the  desire  to 
keep  the  rotational  energy  below  the  energy  required  to  excite  the  first  vibrational  mode 
of  the  H2  molecule. 
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Appendix  F  The  Adiabatic  PESs  of  B  +  H2 


In  this  section  the  adiabatic  PESs  of  the  B  +  H2  system  are  presented  for  the 
configurations  with  C?v  symmetry  (9  =  n/2)  and  Cmil  symmetry  (9  =  0).  The  PESs  are 
presented  as  a  function  of  r  and  R.  Each  surface  is  rendered  in  a  2D  color  plot  which 
captures  the  entire  surface  at-a-glance.  Each  2D  plot  is  accompanied  by  a  3D  rendering 
to  assist  in  visualizing  the  shape  of  the  surface.  The  PESs  were  rendered  on  a  grid  601  x 
601  points  using  cubic  spline  interpolation.  For  energy  values  greater  than  0.2707  au  the 
PESs  display  wavy  features  caused  by  ringing  from  the  cubic  spline  interpolation.  The 
PES  values  were  mapped  to  a  color  palette  containing  200  values  with  bins  ranging  from 
the  minimum  value  of  the  PES  to  the  maximum  value.  The  color  scale  of  the  2D  and  3D 
renderings  are  identical  for  a  given  PES  orientation.  The  figures  begin  on  the  next  page. 
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Figure  46.  The  2D  color  rendering  of  the  adiabatic  PES  =  0) 
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Figure  47.  The  3D  rendering  of  the  adiabatic  PES  =  0) 
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Figure  48.  The  2D  color  rendering  of  the  adiabatic  PES  V-pAi(fi  —  n/2) 
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Figure  49.  The  3D  rendering  of  the  adiabatic  PES  V^zAi(d  =  n/2) 
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Figure  50.  The  2D  color  rendering  of  the  adiabatic  PES  V2zAr(0  —  0) 
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Figure  51.  The  3D  rendering  of  the  adiabatic  PES  V2zAr(0  =  0) 
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Figure  52.  The  2D  color  rendering  of  the  adiabatic  PES  V22A>(9  —  n/2) 
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Figure  53.  The  3D  rendering  of  the  adiabatic  PES  V2-i  A'{9  =  n/2) 
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Figure  54.  The  2D  color  rendering  of  the  adiabatic  PES  VA"(Jd  —  0) 
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Figure  55.  The  3D  rendering  of  the  adiabatic  PES  VA"(d  =  0) 
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Figure  56.  The  2D  color  rendering  of  the  adiabatic  PES  VA"{0  —  tr/2) 
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Figure  57.  The  3D  rendering  of  the  adiabatic  PES  VA"(0  =  n/2) 
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Appendix  G  The  Derivative  Coupling  Surfaces  of  B  +  H2 


In  this  section  the  components  of  the  derivative  coupling  vector  field  are 
displayed  as  a  function  of  r  and  R  for  9  =  n/20.  To  illustrate  the  jagged  features  that 
can  exist  in  DCTs  no  attempt  is  made  to  interpolate  these  surfaces. 


Figure  58.  The  derivative  coupling  surface  for  the  component  12  for  6  —  n/20 


134 


R,  1 2  R,12 


Figure  59.  The  derivative  coupling  surface  for  the  component  r^12  for  6  =  n/20 


Figure  60.  The  derivative  coupling  surface  for  the  component  T<^12  for  8  —  n/20 
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Appendix  H  The  Diabatic  PESs  of  B  +  H2 


In  this  section  the  diabatic  PESs  of  the  B  +  PE  system  are  presented  for  0  =  7r/4. 
The  PESs  are  presented  as  a  function  of  r  and  R.  Due  to  the  symmetry  derived  boundary 
conditions,  diabatic  PESs  for  the  configurations  with  C2V  symmetry  (0  =  n/2)  and  Dooh 
symmetry  (0  =  0)  are  identical  to  their  adiabatic  counterparts.  The  diabatic  coupling 
surface  Vxz  is  zero  for  these  configurations  as  well.  Furthermore,  the  single  surface  with 
A”  Cs  symmetry  does  not  mix  with  other  surfaces.  As  a  result  the  adiabatic  VAn  and 
diabatic  Vyy  PES  are  equal. 

Each  surface  is  rendered  in  a  2D  color  plot  which  captures  the  entire  surface  at-a- 
glance.  Each  2D  plot  is  accompanied  by  a  3D  rendering  to  assist  in  visualizing  the  shape 
of  the  surface.  The  PESs  were  rendered  on  a  grid  601  x  601  points  using  cubic  spline 
interpolation.  For  energy  values  greater  than  0.2707  au  the  PESs  display  wavy  features 
caused  by  ringing  from  the  cubic  spline  interpolation.  The  PES  values  were  mapped  to  a 
color  palette  containing  200  values  with  bins  ranging  from  the  minimum  value  of  the  PES 
to  the  maximum  value.  The  color  scale  of  the  2D  and  3D  renderings  are  identical  for  a 
given  PES  orientation.  The  figures  begin  on  the  next  page. 
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Figure  61.  The  2D  color  rendering  of  the  diabatic  PES  Vxx(8  =  7r/4) 
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Figure  62.  The  3D  rendering  of  the  diabatic  PES  Vx  x(6  —  n/\ ) 
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Figure  63.  The  2D  color  rendering  of  the  diabatic  PES  Vzz(9  =  7r/4) 
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Figure  64.  The  3D  rendering  of  the  diabatic  PES  Vzz{9  —  tc/4) 
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Figure  65.  The  2D  color  rendering  of  the  diabatic  PES  Vxz(8  =  tt/4) 
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Figure  66.  The  3D  rendering  of  the  diabatic  PES  Vxz(8  —  7r/4) 
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Appendix  I  The  Diabatic  Effective  PESs  of  B  +  H2 


In  this  section  a  selection  of  diabatic  effective  PESs  of  the  B  +  H2  system  are 
presented.  The  PESs  are  presented  as  a  function  of  r  and  R.  Each  surface  is  rendered  in 
a  2D  color  plot  which  captures  the  entire  surface  at-a-glance.  Each  2D  plot  is 
accompanied  by  a  3D  rendering  to  assist  in  visualizing  the  shape  of  the  surface.  The 
PESs  were  rendered  on  a  grid  601  x  601  points  using  cubic  spline  interpolation.  For 
energy  values  greater  than  0.2707  au  the  PESs  display  wavy  features  caused  by  ringing 
from  the  cubic  spline  interpolation.  The  PES  values  were  mapped  to  a  color  palette 
containing  200  values  with  bins  ranging  from  the  minimum  value  of  the  PES  to  the 
maximum  value.  The  color  scale  of  the  2D  and  3D  renderings  are  identical  for  a  given 
PES  orientation.  The  figures  begin  on  the  next  page. 
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Figure  67.  The  2D  color  rendering  of  the  diabatic  effective  PES  j  0, 0,  j,  ^  tyeff  |o,  0 
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Figure  68.  The  3D  rendering  of  the  diabatic  effective  PES  (0, 0,  ^  |Pe//  0,  i 
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Figure  69.  The  2D  color  rendering  of  the  diabatic  effective  PES  (4, ~ \Veff  \^>  1, 
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Figure  70.  The  3D  rendering  of  the  diabatic  effective  PES  (4,  |  |  4, 


142 


0.2496 
0.2201 
0.1905 
0.161 
0.1315 
0.102 
0.0725 
0.043 
0.0135 

01  23456789  10 

R  (au) 

Figure  71.  The  2D  color  rendering  of  the  diabatic  effective  PES  ^6,  —  1,  ^  |Pe//|6, 
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Figure  72.  The  3D  rendering  of  the  diabatic  effective  PES  (6, 
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Figure  73.  The  2D  color  rendering  of  the  diabatic  effective  PES 
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Figure  74.  The  3D  rendering  of  the  diabatic  effective  PES  (0, 0,  j,  ^  |Pe//  0,  i 
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Figure  75.  The  2D  color  rendering  of  the  diabatic  effective  PES  (O,  ^,\,\\Veff 
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Appendix  J  The  Adiabatic  Effective  PESs  of  B  +  H2 


In  this  section  a  selection  of  adiabatic  effective  PESs  of  the  B  +  H2  system  are 
presented.  The  PESs  are  presented  as  a  function  of  r  and  R.  Each  surface  is  rendered  in 
a  2D  color  plot  which  captures  the  entire  surface  at-a-glance.  Each  2D  plot  is 
accompanied  by  a  3D  rendering  to  assist  in  visualizing  the  shape  of  the  surface.  The 
PESs  were  rendered  on  a  grid  601  x  601  points  using  cubic  spline  interpolation.  For 
energy  values  greater  than  0.2707  au  the  PESs  display  wavy  features  caused  by  ringing 
from  the  cubic  spline  interpolation.  The  PES  values  were  mapped  to  a  color  palette 
containing  200  values  with  bins  ranging  from  the  minimum  value  of  the  PES  to  the 
maximum  value.  The  color  scale  of  the  2D  and  3D  renderings  are  identical  for  a  given 
PES  orientation.  The  figures  begin  on  the  next  page. 
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Figure  79.  The  2D  color  rendering  of  the  first  adiabatic  effective  PES 
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Figure  80.  The  3D  rendering  of  the  first  adiabatic  effective  PES 
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Figure  81.  The  2D  color  rendering  of  the  10th  adiabatic  effective  PES 


0.2499 

0.2198 

0.1896 

0.1595 

0.1293 

0.0992 

0.069 

0.0389 

0.0087 


Figure  82.  The  3D  rendering  of  the  10th  adiabatic  effective  PES 
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Figure  83.  The  2D  color  rendering  of  the  20th  adiabatic  effective  PES 
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Figure  84.  The  3D  rendering  of  the  20th  adiabatic  effective  PES 
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